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NUMERICAL  AND  ANALYTICAL  STUDIES  OF  STEFAN  PROBLEMS 
FINAL  Technical  Report 


ABSTRACT: 

Compact  finite  element  scheme  is  used  to  solve  Stefan  Problem 
with  one-dimensional  and  two-dimensional  geometries  numerically,  us 
ing  an  enthalpy  formation.  The  numerical  results  indicate  that  the 
position  of  the  melting  front  can  be  determined  with  first  order  ac 
curacy  by  this  method,  the  number  of  iterations  at  each  time  step 
being  determined  largely  by  the  number  of  cells  traversed  by  the 
front  during  a  time  step. 

The  codes  for  both  1 -D  and  2-D  Stefan  problems  are  written  in 
Cray  Fortran  and  vectorization  on  the  Cray  Y-MP  was  used.  The  en¬ 
thalpy  formulation  of  the  1 -D  and  2-D  Stefan  problems  are  approxi¬ 
mated  by  compact  schemes.  The  numerical  results  are  compared  to 
known  exponential  solutions,  and  the  solutions  and  errors  are  plot¬ 
ted  using  mathematica. 

Four  papers  have  been  published  or  completed  for  publication. 
Three  faculty  members,  Bolindra  N.  Borah  (P.I.),  Robert  E.  White 
(Co-P.I.),  and  Milton  E.  Rose  (Co-P.I.)  did  work  in  this  project. 
Besides  three  professors,  there  were  three  graduate  students  who 
also  helped  to  complete  the  project. 
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A  COMPREHENSIVE  LIST  OF  THE  OBJECTIVES  OF  RESEARCH  EFFORT  OR  THE 
STATEMENT  OF  WORK: 


This  research  aimed  at  developing 

(1)  Compact  finite  volume  schemes  for  treating  one  and  two  dimen¬ 
sional  Stefan  Problems  for  bodies  of  general  shapes  and 
connectivities , 

(2)  Analytical  approximations  for  the  position  of  the  interface 
in  problems  with  cylindrical  and  spherical  symmetry,  and 

(3)  Development  and  testing  of  parallel  algorithms  for  nector/ 
multiprocessing  computers. 
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2.  STATUS  OF  THE  RESEARCH  EFFORT:  A  substantive  statement  of  signifi¬ 
cant  accomplishments  and  progress  towards  achieving  the  research 
objectives 

The  compact  finite  volume  method  was  used  to  treat  general 
diffusion  equations.  By  using  intrinsic  geometrical  properties  of 
the  volume  element,  it  was  possible  to  describe  discrete  versions 
of  the  div,  curl  and  grad  operations  which  lead,  using  summation- 
by-parts  techniques,  to  familiar  energy  equations  as  well  as 
divcurl  =  0  and  curldiv  =  0  identities.  For  the  diffusion  equa¬ 
tions,  these  operators  describe  compact  schemes  whose  convergence 
is  assured  by  the  energy  equations  and  which  yield  both  potential 
and  flux  vectors  with  second  order  accuracy.  A  simplified  poten¬ 
tial  form  is  specially  useful  for  obtaining  numerical  results  by 
multigrid  and  ADI  methods.  The  treatment  of  general  curvilinear 
coordinates  was  shown  to  result  from  a  application  of  these  gen¬ 
eral  results. 

An  implicit  enthalpy  scheme  is  used  to  solve  one-dimensional 
and  2-D  Stefan  problems  by  compact  finite  volume  scheme.  The  re¬ 
sults  are  compared  with  known  analytical  solutions  (exponential 
solution)  and  the  solutions  and  the  errors  are  plotted  using  Mathe- 
matica.  Many  intresting  graphs  are  presented.  Plot  of  the  numeri¬ 
cally  calculated  front  positions  and  the  least  square  fit  for  these 
numerically  determined  front  positions  are  supplied.  In  two  dimen¬ 
sional  problems  the  3-D  plot  of  the  numerical  solution  shows  the 
behavior  in  the  neighborhood  of  the  front.  A  3-D  plot  of  the  numeri¬ 
cal  error  shows  that  the  larger  errors  are  mostly  confined  to  the 
neighborhood  of  the  front. 
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The  2-D  Stefan  Problem  is  solved  for  both  Dirichlet  and  Neumann 
boundary  conditions. 

The  parallel  algorithms  for  the  enthalpy  formulation  of  multi¬ 
ple  dimensional  two  phase  Stefan  problem  have  been  studied.  The 
compact  scheme  was  shown  to  be  adaptable  to  standard  parallel  al¬ 
gorithmic  techniques  for  the  efficient  vectorization  of  multidimen¬ 
sional  Stefan  Problems. 

The  compact-ADI  scheme  for  2-D  Stefan  Problem  was  compared 
with  compact  SOR  and  traditional  SOR  method.  Among  all  three 
method,  compact-ADI  scheme  was  found  to  be  the  most  efficient 
method  as  regards  to  time  of  computation  and  memory  usages.  All 
three  methods  were  computed  in  Alliant  FX/40  and  Cray  Y-MP  vector/ 
multiprocessing  computer.  These  results  were  scheduled  to  be  pub¬ 
lished  in  paper  (IV)  listed  in  section  3. 


3.  A  CUMULATIVE  CHRONOLOGICAL  LIST  OF  WRITTEN  PUBLICATIONS  IN  TECHNI¬ 
CAL  JOURNALS: 

The  following  four  papers  are  either  published  or  ready  for 
publication: 

(i)  Compact  Finite  Volume  Methods  for  Diffusion  Equations,  Milton 
E.  Rose,  Journal  of  Scientific  Computing,  Vol.  4,  No.  3,  pages 
261-290,  Sept.,  1989  (reprint  enclosed). 

(ii)  An  Implicit  Enthalpy  Scheme  for  One-phase  Stefan  Problems, 
Milton  E.  Rose,  submitted  for  publication  in  the  Journal 
of  Scientific  Computing  (manuscript  enclosed). 

(iii)  A  Comparative  Studies  of  Compact  Volume  Methods  for  the  2-D 
Diffusion  Equation  with  Finite  Difference  ADI  and  SOR. 

Bolindra  N.  Borah,  Robert  E.  White  and  Archimedes  Kyrillidis. 

(iv)  A  Compact  Finite  Volume  Scheme  for  2-D  Stefan  Problems  and 
Vector/Multiprocessing  Computers,  Robert  E.  White,  Bolindra 
N.  Borah  and  Archimedes  Kyrillidis  (ready  for  publication, 
manuscript  enclosed). 
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4.  LIST  OF  PROFESSIONAL  PERSONNEL  ASSOCIATED  WITH  THE  RESEARCH: 

The  following  faculty  members  have  worked  in  this  project: 

(i)  Bolindra  N.  Borah,  Professor  of  Applied  Mathematics  and  Com¬ 
puter  Science  (P.I.),  North  Carolina  A&T  State  University, 
Greensboro,  NC. 

(ii)  Milton  E.  Rose,  Adjunct  Professor  of  Mathematics  and  Com¬ 
puter  Science  (Co-P.I.),  North  Carolina  A&T  State  University, 
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(iii)  Robert  E.  White  (Consultant)  Professor  of  Mathematics,  North 
Carolina  State  University,  Raleigh,  NC. 

The  following  graduate  students  have  worked  in  this  projects: 
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5.  INTERACTIONS:  PAPER  PRESENTED  AT  MEETINGS,  CONFERENCES,  SEMINARS, 
ETC. 

(i)  The  paper  entitled  "A  Compact  ADI  Scheme  for  2-D  Stefan  Pro¬ 
blems  and  Vector/Multiprocessing  Computers"  will  be  presented 
at  the  Second  International  Conference  on  Industrial  and  Ap¬ 
plied  Mathematics,  July  8-12,  1991  to  be  held  at  Washington, 

D.  C. 

(ii)  The  entire  research  findings  were  presented  at  the  seminar 

organized  by  the  McDonnell  Douglas  Aircraft  Company,  St.  Louis, 
MO,  63166  held  April  23,  24,  1991. 
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Compact  Finite  Volume  Methods 
for  the  Diffusion  Equation 

Milton  E.  Rose‘ 

Received  September  21.  1989 


We  describe  an  approach  to  treating  initial-boundary-value  problems  by  finite 
volume  methods  in  which  the  parallel  between  dilTerential  and  difference 
arguments  is  closely  maintained.  By  using  intrinsic  geometrical  properties  of  the 
volume  elements,  we  are  able  to  describe  discrete  versions  of  the  div.  curl,  and 
grad  operators  which  lead,  using  summation-by-parts  techniques,  to  familiar 
energy  equations  as  well  as  the  div  curl  =  0  and  curl  grad  =  0  identities.  For  the 
diffusion  equation,  these  operators  describe  compact  schemes  whose  con¬ 
vergence  is  assured  by  the  energy  equations  and  which  yield  both  the  potential 
and  the  flux  vector  with  second-order  accuracy.  A  simplified  potential  form  is 
especially  useful  for  obtaining  numerical  results  by  multigrid  and  ADI  methods. 
The  treatment  of  general  curvilinear  coordinates  is  shown  to  result  from  a 
specialization  of  these  general  results. 


KEY  WORDS:  Finite  volumes;  compact  schemes;  elliptic;  diffusion  equation; 
curvilinear  coordinates. 

1.  INTRODUCTION 

Let  L  be  a  domain  with  boundary  surface  S  on  which  n  is  the  outward  unit 
normal.  This  article  describes  a  finite  volume  scheme  for  solving  the 
diffusion  equation  for  a  potential  <t>  and  flux  u  in  the  form 


D.E. 

=  div  u  —  /  in  L 

(1.1) 

grad  ^  =  u 

B.C. 

otvi  ■  n  +  fltf)  =  g  on  5,  /  ^  0 

(1.2) 

I.C. 

<j>  =  go  in  E,  /  =  0 

(1.3) 

'  Department  of  Mechanical  Engineering,  N.C.  A&T  State  University,  Greensboro,  North 
Carolina  27411.  Present  Address:  4505  Tower  Road,  Greensboro,  North  Carolina  27410. 
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in  both  the  steady  and  unsteady  cases.  We  recall  that  the  solution  of  this 
problem  satisfies  an  “energy”  equation 


d_ 

7t 


(I>\IV  + 

V  r 


u  •  u  dV  = 


(1.4) 


which  follows  by  multiplying  the  first  equation  in  (1.1 )  by  and  integrating 
by  parts.  We  also  recall  that  when  f  =0  the  maximum-minimum  values  of 
(j)  either  lie  on  the  boundary  S  or,  in  case  /  =  0,  in  V  itself. 

A  discrete  form  of  these  results  will  hold  for  the  finite  volume  scheme 
and  includes  results  about  generalized  curvilinear  coordinate  mappings 
as  a  special  case.  We  will  identify  finite  volume  operators  div;,,  grad;,, 
and  curl,,  which  are  consistent  with  their  differential  counterparts  and 
from  which  a  discrete  energy  expression  corresponding  to  (1.4)  will 
follow  by  simple  summation-by-parts  identities.  The  important  identities 
div,,  curl,,  =  0  and  curl,,  grad,,  =  0  will  also  remain  valid,  as  we  will  show  at 
a  later  point  in  the  article.  A  potential  form  involving  the  operator 
div,,  grad,,  which  is  obtained  by  eliminating  u,  leads  to  a  symmetric, 
positive  definite  operator  to  which  multigrid  and  other  fast  solution  techni¬ 
ques  are  applicable. 

The  fact  that  the  finite  volume  schemes  described  here  lead  to  discrete 
energy  expressions  is  the  principal  result  of  this  article.  It  ensures  that  the 
schemes  converge.  We  will  also  find  that  the  approximations  to  both  (p  and 
u  will  be  second  order  accurate.  This  result  is  similar  to  that  obtained  for 
mixed  finite  element  methods  (Brezzi,  1988;  Strang  and  Fix.  1973)  and  it 
is  possible,  in  fact,  to  view  the  present  scheme  as  a  finite  element  method 
that  involves  nonconforming  elements  (Gatski  et  ai,  1989;  Lustman  and 
Rose,  1988). 

Many  of  these  ideas  can  be  illustrated  most  simply  for  steady,  one¬ 
dimensional  problems,  for  which  reason  we  first  discuss  the  equations 
=  u,  u'  =  /  in  detail.  We  will  find  it  natural  to  introduce  a  primary  mesh, 
which  is  formed  by  the  end  points  of  subintervals  into  which  the  basic 
domain  is  divided,  and  a  dual  mesh,  which  is  formed  by  the  center  points 
of  the  primary  mesh.  The  variables  p  and  u  will  be  associated  with  the 
primary  mesh  while  another  variable  {(/  will  be  associated  with  the  dual 
mesh.  An  algebraic  relationship  between  p,  0,  and  u  on  each  subinterval 
provides  an  approximation  to  the  solution  operator  (for  which  reason  the 
scheme  is  called  compact)  and  the  solution  in  the  large  is  obtained  by 
imposing  continuity  conditions  for  <P  and  u  at  points  of  the  primary  mesh. 
Both  <f>  and  i}/  will  be  found  to  converge  with  second-order  accuracy  to  the 
solution  at  the  points  at  which  they  are  defined.  Furthermore,  as  noted 
above,  although  u  will  be  defined  by  one-sided  divided  differences  involving 
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(f>  and  ip,  its  convergence  will  also  be  second  order  accurate  as  a  result  of 
the  continuity  conditions  imposed. 

In  extending  these  ideas  to  higher-dimensional  problems  by  sub¬ 
dividing  the  solution  domain  V  into  volume  elements  we  will  also  find  that 
a  primary  and  dual  grid  play  a  natural  role.  Now  the  variables  (f>  and  u  will 
be  associated  with  the  center  points  of  the  faces  of  a  volume  element  while 
if/  will  be  associated  with  the  center  of  the  element.  However,  an  additional 
variable  (  (or  a  box-variable  y)  which  is  associated  with  the  edges  (ver¬ 
tices)  of  the  element  will  also  be  required.  The  compact  scheme  will 
describe  relationships  between  these  variables  which  produce  second-order 
accurate  results  at  the  points  indicated.  In  the  method  discussed  here,  the 
nontangential  components  of  u  on  a  face  of  an  element  are  obtained  with 
the  use  of  one-sided  differences  in  terms  of  (P  and  ip  while  its  tangential 
components  are  determined  solely  by  the  edge  values  C.  and  these  are 
obtained  in  terms  of  ip  at  points  of  the  dual  grid  by  bilinear  interpolation. 
In  the  case  of  (uniform  or  nonuniform)  Cartesian  or  orthogonal  coordinate 
elements,  the  variables  C  as  well  as  the  tangential  components  of  u  on  faces 
of  elements  can  be  obtained  by  postprocessing,  if  required. 

Some  of  these  ideas  are  familiar  from  applications  of  finite  volume 
methods  to  fluid  dynamics  and  are  described  in  Peyret  and  Taylor  (1983) 
and  in  a  review  paper  by  Vinokur  (1989).  A  focus  of  many  such  methods 
is  on  the  treatment  of  conservation  laws  for  fluids.  As  noted  by  these 
authors,  the  primary  and  dual  grids  often  play  a  role  in  many  such  schemes 
whenever  gradient  terms  are  included,  as  occur  when  viscous  effects  arise. 
Indeed,  the  need  to  accurately  approximate  quantities  like  the  stress  tensor 
at  boundaries  of  general  domains  is  a  principal  reason  to  resort  to  finite 
volume  methods,  although  their  use  may  add  significantly  to  the  cost  of 
computations.  To  understand  some  of  the  problems  that  can  arise,  the 
diffusion  equation  can  serve  as  a  useful  model.  We  will  find  that  the 
relationships  between  grid  variables  dificr  in  important  respects  from  those 
described  elsewhere.  Many  schemes  place  primary  emphasize  on  the  vertex 
(box)  variables  and  most  methods  deliberately  avoid  the  use  of  one-sided 
differencing  to  approximate  the  flux  except,  perhaps,  at  boundaries. 
Although  compact  schemes  related  to  those  described  here  have  been  used 
for  Cartesian  grids,  e.g.,  Gatski  ei  al.  (1989)  and  Rose  (1983),  the  roles  of 
the  variables  were  not  completely  developed. 

The  paper  concludes  by  describing  the  rather  straightforward 
modifications  that  arc  required  to  treat  the  time-dependent  problem.  The 
result  is  a  Crank  Nicholson-type  scheme  and  energy  arguments  provide 
convergence  estimates  for  the  finite  volume  method.  We  arc  also  able  to 
show  that  energy  arguments  can  be  adapted  to  a  Peaceman  -Rachford  ADI 
scheme. 
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2.  STEADY  ONE-DIIMENSIONAL  PROBLEMS 

2.1.  Notations 

On  an  interval  [I  ,  I  ^  ]  consider  the  problem; 

D.E.  »'  =  /  (2.1a) 

(j)'  =  u 

B.C  </>=g  for  .v=l,  (2.1b) 

The  associated  energy  equation  is 

(f)f  (lx  +  j  ir  (lx  =  {(f>u)^^-{(l>u)^  ( 2.2 ) 

Divide  [1  _  ,  1  ^.  ]  into  M  nonoverlapping  subintervals  =  {-v|-v^  _  ,  2  ^ 
-'c  \/2}  center  points  .Vy,  y  =  1/2,  3/2,...,  M  —  1/2.  The  points  .v,, 

7=1,  2,...,  M  —  1,  are  end  points  of  the  subintervals  which  lie  interior  to 
[1  ,  1  +  ]-  We  let  /,,,  4  denote  the  sets  of  indices  corresponding  to  these 
points; 

/,.=  {!,  2,...,  M  -  1 }  (interior  end  points) 

4  =  {1/2,  3/2,...,  A/  -  1/2 }  (center  points ) 

We  adopt  the  finite-difference  notations  ^.v,  =  .y,+ ,^2  — -v,  1  2-  /?,  =  ‘d.v,/2, 
u(x,)  =  Ui-  By  introducing  the  central  average  and  difference  operators 

^L(t>i={<l>i+M2  +  (f>i-\nMX  Mi={<l>i^\:2-<t>,-\2)  (2.3) 

we  can  verify  the  summation-by-parts  identities 

A{(l>\lj)  =  {^i((>){A4/)  + {A(t>){^i\lf)  (2.4a) 

/<(^i/^)  =  (/r^)(/4)  +  (/l<!))(-di/^)/4  (2.4b) 

Both  will  play  a  central  role  in  establishing  energy  results. 

2.2.  A  Compact  Scheme 

In  each  subinterval  /,,  construct  an  approximate  solution  using  values 
,  1/2  boundary  data.  Specifically,  (i)  for  Je  /,  set 


(2.5) 

flA*,  1.2- ^j) 

(2.6a) 

h,U,  =  |.) 

(2.6b) 
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4),  U  (j).  U 

Fig.  1.  An  interval  showing  the  association  of  the  variables  u  with  its  end  points  and  of 

i//  with  its  center  point. 


Adopting  the  convention  that  indicates  an  approximation  to  the  potential 
solution  at  the  center  of  the  interval  while  +  1/2  indicates  the  approximation 
as  its  end  points  we  see  that  (2.5)  is  a  central  difference  approximation  to 
\x'—f  while  (2.6)  approximates  =  u  by  one-sided  differences.  Then 
Eqs.  (2.5)  and  (2.6)  can  be  solved  for  ^2  and  \}/j  in  each  interval  in  terms 
of  (l>j+  1/2  considered  as  boundary  data. 

(ii)  Next,  require  that  both  u  and  tp  shall  be  continuous  across  every 
end  point  common  to  two  intervals,  i.e., 

[i/],  =  [^],  =  0  for  /e  4  (end  points)  (2.7) 

(This  is  implied  by  our  notation,  since  we  have  not  distinguished  between 
the  right-  and  left-hand  limits  of  u  and  <f>  at  end  points.)  Using  (2.6)  to 
evaluate  m,+  i/2  at  the  right  and  left  end  points  of  the  adjacent  intervals 
4-  1/2,  4+  1/2  we  find 

^-v,  +  i/2W,+  ,/2  =  ('A,+  i-i/',-iK  'e4>  (2.8) 

which  indicates  that  the  accuracy  of  u  may  be  higher  than  that  suggested 
by  the  one-sided  difference  expressions  which  originally  defined  u  in  (2.6). 
Also  note  that  (2.6)  can  be  used  to  express  the  continuity  condition 
[«],  =  0  in  terms  of  <f>,  ip  with  the  result 

((pj  —  ip , _  \  i2)lhj  1/2  =  (•A/  + 1/2  “  +  1/2  (2.9a ) 

and  which,  in  the  case  of  equal  intervals,  reduces  to 

(pi  =  fup,,  /e /.  (end  points)  (2.9b) 

The  addition  of  the  boundary  conditions  leads  to  a  determined  system  of 
algebraic  equations  for  the  M  values  of  ip  and  the  M  -f-  1  values  of  (p  and  u. 

The  special  nature  of  this  system  is  worth  noting.  By  eliminating  u 
from  (2.5)  using  (2.6)  we  find,  in  the  case  of  equal  intervals, 

{p<p ,  -  ip ,)  =  2hf  \  /G /,  (center  points)  (2.10a) 
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and  if  we  also  write  the  continuity  conditions  (2.9b)  in  the  form 

=  0,  /e /,.  (end  points)  (2.10b) 

it  is  evident  that  t//,  are  the  odd-even  components  of  the  solution  vector 
of  a  tridiagonal  system  of  equations  which  is,  thus,  easily  solved.  In  sharp 
contrast  to  (2.10),  the  standard  finite  difference  scheme  for  solving  (2.1 )  is 

/e/.  u/,  (2.10c) 

so  that  some  differences  with  standard  numerical  treatments  of  differential 
equations  are  to  be  expected  using  the  schemes  considered  in  this  article. 


2.3.  The  Energy  Equation  and  Convergence 

We  will  now  show  that  this  scheme  leads  to  a  discrete  form  of  the 
energy  equation  (1.4).  For  simplicity,  we  will  assume  a  uniform  m  sh.  Add 
and  subtract  the  expressions 

{Axil)  =  ±((;^^±  ,,,  -  I//,) 

given  in  (2.6)  to  obtain 

Axnu,  =  A<f)^  (2.11a) 

ij/j  =  — Ax  AUi/4  (2.11b) 

and  we  see  that  the  approximations  to  ^  and  »//  are  second  order  with 
respect  to  the  center  of  an  interval.  Then,  starting  with 

AUi^Axfi  (2.5) 

multiply  by  ip  and  use  (2.11)  and  (2.4)  to  obtain 

^  Au  =  {pip  —  {Ax/4)  ziu)  ■  An 

=  A{(l>u)  —  {pu  A<p  +  (d.Y/4 )(.^h/)" ) 

=  A{(Pu)  —  /lv[  (/<!/)"  +  (/l//)‘/4  J 
=  A{<f>u)  —  Ax  p{ir)  (2.12) 


(Here  and  in  the  following  we  omit  subscripts  when  no  confusion  is  likely 
to  arise.)  Summing  over  center  points  of  intervals  we  obtain 

1 1 . 


X  (<A,  AUi  +  Ax  pu^ )  =  {<l>u) 


/, 


/, 
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Because  of  the  continuity  conditions,  the  jumps  [  ]  vanish  and  we  obtain 
the  discrete  energy  expression 

I . 

X  {'A/.//  +  /'"/■)  ^  (2.1 3) 

/,  I 

'('his  result  will  enable  us  to  conclude  that  (ft,  xfj,  and  u  converge  with 
second-order  accuracy.  Introduce  center-point  and  end-point  norms  as 
follows; 

=  (2-l4a) 

lii'll,.=  +  +  (2.l4b| 

By  summing  (2.8)  and  using  (2.6)  we  find 

-^0=  H  +  X  *6 

A- 

and  if  we  also  assume  end-point  conditions  <f>o  =  </>ju  =  ^  we  obtain  the 
inequality 

The  energy  expression  itself  leads  to  the  estimate 

ll/ll,. 


while  the  continuity  condition  (2.9)  results  in 


Thus,  we  have 

ll/ll, 

i!u||,,^C||/|i,,  (2.15) 

ll/ll, 

Applying  this  to  the  homogeneous  problem  we  conclude  that  u  =  (ft  ~ 
iff  ={).  Thus  the  algebraic  problem  has  a  solution  that  is  unique.  Next, 
interpreting  /  as  the  truncation  error,  we  have /=  f7(/1.v“),  and  we  see  that 
(ft,  if/,  and  II  converge  with  second-order  accuracy. 

This  conclusion  about  the  second-order  accuracy  of  u  applies  not  only 
to  its  values  at  interior  end  points,  as  might  be  expected  from  (2.8),  but 
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Table  I.  Error  Norms 


Number  of 
intervals 

6 

l.234e-2 

12 

3.472e-3 

9.259e-3 

24 

8.680<'-4 

2.893e-3 

48 

2.l70t>-4 

7.957e-4 

96 

5.425t'-5 

2.()S0('-4 

192 

l.356e-5 

5.312e-5 

also  to  the  values  Uq  and  u  ,^,  at  the  end  points  of  [1  .  1  ]  itself  and 

which,  we  recall,  were  defined  in  Eq.  (2.6)  by  one-sided  difference  expres¬ 
sions. 

Table  I  presents  computations  which  verify  these  conclusions  for  the 
differential  equation  (2.1)  corresponding  to  the  solution  0  =  .v‘(l— .v)^ 
Note  that  the  errors  are  measured  at  the  end  points  of  intervals. 

Finally,  we  remark  that  more  general  boundary  conditions  involving 
both  u  and  <f>  also  lead  to  the  same  conclusions. 


2.4.  The  Potential  Form 

By  eliminating  the  flux  u  in  (2.1)  we  obtain  the  familiar  second-order 
equation  (f>"  =  /which  we  call  the  potential  form.  The  difference  scheme  that 
corresponds  to  this  may  be  obtained  by  eliminating  u  in  the  scheme  just 
described.  To  illustrate,  using  (2.6),  eliminate  u  in 

AUj  =  /i.\-J]  (2.5) 


to  obtain 


h  =  Ax  12 


(2.16) 


As  we  have  seen,  for  uniform  meshes  the  continuity  conditions  [u]=0 
result  in 


=  iel,.  (2.9b) 

To  these  are  to  be  added  the  boundary  conditions 


(/>ii  =  ^  > 
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When  /  =0  we  obtain 

^  =  (2.17) 

and  the  maximum-minimum  property  of  the  solution  is  an  immediate  con¬ 
sequence.  This  can  be  used  to  develop  an  alternate  proof  of  convergence. 

The  potential  form  for  the  discrete  problem  can  be  seen  to  arise  as  a 
variational  problem  associated  with  the  discrete  energy  equation.  The 
maximum-minimum  property  reflects  the  fact  that  the  associated  algebraic 
problem  is  symmetric  and  positive  definite.  In  the  general  treatment,  the 
potential  form  will  be  more  suitable  for  numerical  work,  while  the  general 
form  involving  <}>,  if/,  and  w  as  a  first-order  system  will  prove  more  con¬ 
venient  for  theoretical  discussions,  particularly  for  the  development  of 
energy  estimates. 

3.  A  COMPACT  FINITE-VOLUME  METHOD 

We  now  turn  to  the  problem  of  treating  the  boundary  value  problem 
corresponding  to  the  steady-state  solution  of  (1.1)  when  F  is  a  general 
volume.  Our  objective  will  be  to  partition  V  into  volume  elements  <5F  in 
such  a  manner  that  the  prescribed  boundary  data  on  S  can  be  accurately 
transmitted  to  the  boundary  elements  and  then,  by  solving  a  discrete 
boundary  value  problem  in  each  element  corresponding  to  a  compact 
scheme,  obtain  an  approximate  solution  in  all  elements  which  also  satisfies 
an  energy  expression,  thus  ensuring  convergence.  As  a  result,  we  expect  to 
be  able  to  treat  problems  posed  in  curvilinear  coordinates  as  a  special,  but 
important,  case. 

In  order  to  be  able  to  generalize  the  arguments  given  earlier,  a  number 
of  additional  notations  will  have  to  be  introduced.  Before  doing  so,  we  will 
first  present  a  short  overview  of  the  principal  ideas  that  will  be  involved. 
The  more  detailed  discussion  and  demonstration  that  the  final  result  again 
leads  to  an  energy  expression  and  thus  produces  a  convergent  scheme  may 
be  omitted  if  the  reader  wishes  to  turn  to  the  discussion  of  the  time- 
dependent  problem  given  in  Section  4. 

3.1.  An  Overview 

Guided  by  the  earlier  discussion  of  the  one-dimensional  problem,  we 
may  identify  the  following  requirements  to  construct  a  discrete  approxima¬ 
tion  to  div  u  =  /,  u  =  grad  <f>  when  volume  elements  are  involved: 

i.  Construct  consistent,  discrete  approximations  to  div  u,  grad  (}>  in 
a  volume  element  SV  in  terms  of  variables  associated  with  u  and 
<f/  at  appropriate  points  of  SV. 
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ii.  Further  relate  these  variables  as  may  be  required  so  as  to  lead  to 
a  determined  and  consistent  system  of  algebraic  equations  from 
which  u  and  (f)  can  be  determined  whenever  boundary  conditions 
arc  prescribed  on  S. 

iii.  Develop  a  difference  calculus  that  allows  summation-by-parts 
identities  to  be  applied  to  volume  elements  and  thereby  lead  to  a 
discrete  energy  expression  from  which  a  convergence  argument 
can  be  constructed. 

We  will  confine  our  attention  to  hexahedral  volume  elements  c^V, 
although  we  will  allow  certain  of  its  vertices  to  coalesce,  thereby  including 
tetrahedra  and  related  elements  in  our  discussion.  Figure  2  indicates  the 
center  point  P  of  SV,  the  center  Q  of  an  oriented  face  <5S,  and  a  repre¬ 
sentative  center  point  R  of  an  edge  of  on  which  T  is  a  vertex.  In  addi¬ 
tion  to  the  variables  if/iP)  introduced  in  our  discussion  of  the  one¬ 

dimensional  problem,  we  will  also  associate  another  variable  C{R)  with 
points  R.  These  variables  will  approximate  the  potential  solution  at  the 
points  indicated.  Later,  we  will  also  associate  a  variable  y  Ihe  vertex 
points  T.  Section  3.2  introduces  notations  that  allow  us  to  describe  various 
geometrical  properties  of  an  element  and  from  which  averaging  and 
difference  operators,  corresponding  to  (2.4),  can  be  introduced. 

Section  3.3  describes  discrete  approximations  to  the  operators  div  and 
grad.  The  approximation  to  div  u  is  based  on  the  use  of  Gauss’  theorem  in 
a  volume  element.  Center-point  quadrature  approximations  to  the  surface 
integrals  involved  leads  to  a  discrete  operator  div;,  which  provides  an 
intrinsically  consistent,  second-order  accurate,  approximation  to  div  u  in 
terms  of  the  normal  components  of  u  at  center  points  of  faces  dS.  This  is 


T 


Fig.  2.  A  normal  volume  element;  /’,  Q,  R,  are  the  centers  of  the  volume  i'll',  a  face  (fS.  and 
an  edge  on  which  T  is  a  vertex.  Shown  at  Q  is  a  basis  [<-]  and  <  ihe  unit  normal  to  (>S. 
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a  familiar  idea  in  finite  volume  methods.  An  approximation  to  u  =  grad^ 
will  result  from  the  integral  form 

f  \x  ■  (i\  —  <j){A) 

in  which  the  points  A  and  B  arc  identified  with  the  points  P,  Q.  and  R  in 
Fig.  2.  A  midpoint  evaluation  of  the  integral  will  determine  two  tangential 
components  of  u  lying  in  a  surface  element  in  terms  of  the  variable  C  at 
edges  of  SS{Q).  The  remaining  nontangcntial  component  of  u  on  tlS  will 
result  by  using  points  P  and  Q  and  a  one-sided  evaluation  of  the  integral 
at  Q.  This  approximation  will  lead  to  second-order  accuracy  as  a  result  of 
imposing  additional  continuity  conditions  for  the  variables  on  ^S,  as 
occurred  in  the  one-dimensional  case.  We  indicate  by  the  vector  u 

evaluated  at  each  face  using  this  construction. 

The  compact  scheme  which  emerges  is  described  in  Section  2.4  and 
can  be  summarized  as  follows; 

a.  In  each  element  u,  (f>,  ij/,  (  satisfy  div^  »=/,  grad;,  ^  =  u.  in  which 
the  normal  components  of  u  arising  in  div;,  u  are  to  be  evaluated 
in  terms  of  the  components  given  by  grad;,  <f). 

b.  Across  each  oriented  face  impose  the  continuity  conditions 

in-ds-\  =  w  =  m=o. 

c.  If  3S  lies  on  the  boundary  S  of  F,  then  <j)  and  C  ‘tre  to  be 
prescribed  by  the  data. 

These  conditions  exactly  mirror  the  situation  discussed  for  the  one¬ 
dimensional  problem.  However,  we  shall  find  that  they  lead  to  an  under¬ 
determined  system  of  algebraic  equations  and  certain  additional  conditions 
will  be  required.  These  are  furnished  by  the  completeness  conditions: 

d.  Completeness  Conditions:  These  conditions  determine  c  at  vertex 
points  in  terms  of  the  variable  \J/  by  the  use  of  a  bilinear  interpola¬ 
tion  which  we  indicate  as  l^  =  Q{ij/).  Among  other  things,  this  will 
ensure  that  is  C  bounded  in  norm  by  if/. 

This  scheme  will  lead  to  an  energy  equation  and  convergence  will  be 
assumed.  The  potential  form  will  also  show  that  the  maximum-minimum 
property  also  holds. 

When  curvilinear  coordinates  are  used,  many  of  the  geometrical  con¬ 
structions  required  for  elements  will  be  provided  analytically.  Furthermore, 
for  Cartesian  or  more  general  orthogonal  grids,  wc  will  find  that  the 
variable  C  can  be  obtained  by  postprocessing  the  results  using  the  com¬ 
pleteness  conditions.  This  emphasizes  the  practical  advantages  of  utilizing 
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regular  cleinciils  throughout  nu)st  of  the  dotiiani  I  and  restricting  the 
general  eonstruetion  to  boundary  elements. 

3.2.  Ccomctrical  Considerations  and  isolations 

3.2.1.  Normal  Volume  iUvmcnts 

We  call  jbCj  a  normal  covering  of  V  by  volume  elements  bl',  them¬ 
selves  called  tiormal,  if  (i)each  element  is  a  he.vahedron.  some  of  whose 
face  areas  may  be  allowed  to  vanish  providing  its  volume  i  remains 
positive  {rk'i^cncracy),  (ii)  any  face  of  a  boundary  element  dV  that  is  in 
contact  with  the  boundary  surface  S  of  V  is  tangent  to  S  at  the  center  point 
of  the  face.  The  first  assumption  has  important  consequences  for  the 
geometrical  constructions,  which  will  now  be  described. 

Figure  2  indicates  a  normal  volume  element  with  center  at  P  and 
surface  elements  ^'iS.  Since  SV  is  a  hexahedron,  P  may  be  found  as  the 
average  of  th*'  vertices  of  SV,  while  the  average  of  the  vertices  on  a  face 
yields  the  center  point  Q  and  the  average  of  the  vertex  points  on  an  edge 
yields  the  center  point  R.  Wc  may  thus  construct  a  right-handed,  local 
coordinate  basis 


[e(P)]  =  (e,(TKe,(T).edTn 

at  P  which  is  formed  by  unit  vectors  in  directions  which  connect  the  pairs 
of  opposite  center  points  of  the  faces  OS  which  we  denote  as  0,^, 
i—  1,2,3.  The  face  d^{Q,  )  is  assumed  \>ricnted  into  (')f,  while  dS{Q,,  ) 
is  outward;  these  are  sometimes  called  inHow  and  outflow  faces.  We 
also  assume  that  the  orientation  carries  over  to  neighboring  elements  by 
requiring  that  an  outflow  face  from  one  element  be  an  inflow  face  on  its 
adjacent  neighbor. 

Our  primary  requirement  will  be  to  evaluate  the  vector  u((2)  at  the 
center  points  of  faces.  In  order  to  do  this  we  will  construct  a  unit  basis 
[e((2)J  at  each  Q  as  follows;  assume  that  hS((2)  is  the  outflow  face  of 
dV{P)  and  the  inflow  face  of  a  neighboring  element  c)F(T'),  e.g.,  Q^Q\. 
We  take  e|((2)  as  a  unit  vector  in  the  direction  of  a  line  connectirig  P  and 
P'.  Also  take  eq(_2),  e,((>)  as  unit  vectors  in  directions  that  connect  O  vvith 
the  center  pome  of  the  edges  R  of  <'>S((2)  so  that  [e((2)]  forms  a  right- 
handed  basis  as  indicated  m  the  l  ig.  2.  Thus,  eq(2).  eq^)  are  tangent 
vectors  on  ().S((7|  while  e,((/)  is  nontangential. 

We  may  regard  the  points  indicated  as  being  related  by  displacement 
operators  ,  along  a  direction  c,  so  that 


Q, .  = 
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with  the  further  understanding  tlial  (/,  /,  A  )  indicates  an  even  permutation 
of  (1,  2,  3).  The  averaging  and  difference  operators  denned  earlier  by  (2.3) 
may  be  generalized  by  writing 

=  r)  +  «6(/:,  01/2.  o] 

(3.1) 

for  a  representative  center  point  C.  From  this  it  is  apparent  that  the 
summation-by-parts  identities  (2.4)  continue  to  hold. 

Motivated  by  this  notation  we  will  also  to  indicate  the 

distance  between  the  points  Q,  and  and  write  h,(Q)  —  A ,.\(Q)/2. 
Also,  h  itself  will  indicate  the  minirp.'.'m  value  of  h,(P)  for  all  volume 
elements  covering  V. 

3.2.2.  Covaiioni  and  Conti avarkmi  Ba.se.s  on  SV 

A  basis  [c]  =  (€,,02,03)  may  be  considered  as  a  covariant  basis  from 
which  a  reciprocal  contravariant  basis  [o]  ' '  =  (o',  e\  o^ )  can  be  construc¬ 
ted  satisfying 

e'-(e,  xej/>yc,  o,  =  •  (eO  e^  ),  =  (3.2) 

in  which,  if  c,y  =  o,  •  oond  o'^^o'-o',  then 

e  =  dct(c„),  l/e  =  det(e'^)  (3.3) 

Anticipating  the  connection  to  curvilinear  coordinates,  we  call  element 
metric  coefficient.s  at  P.  As  a  result  we  may  write 

u'  =  ZeX’ 

U  •  II  =  ^  u,u‘ 

where  w,  =  u  o,  and  u'  =  u  c'  are  the  (physical)  covariant  and  contravariant 
components  of  u,  respectively,  with  respect  to  [o]. 

Using  the  unit  tangent  vectors  o^.o^  the  oriented  area  of  a  face 
SSiQ,  ,  )  is 

()S  =  IbSI  •  Oi  X  C3  =  IhS]  ^0  ■  o'  ( 3.5 ) 


where 


IhSI  ^A:.\-iQ)A,.x{Q) 
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Thiis,  the  covarianl  componciil  <>S,  of  <)S  is  dS,  =  Ic^SI  yjc.  Also,  elementary 
geomelrieal  arguments  based  on  the  fact  that  normal  volume  elements  are 
hexahedra  show  that  if 

r„  /l|  A  ■  .1  ..V  .  1 ,  V 

tlien  the  volume  measure  r  of  dl'  is  given  by 

r  =  /t,A-((;)c,(/^)- [<>S((>,  )  +  bS(C>,,  )]  2-r„  //,ye 

i—  I,  2,  3.  If  we  also  assume  positive  constants  c,,  and  c,  such  that  c\^li  ^ 
and  also  set  v^e„  =  //,  ^c  we  can  summarii'c  these  geometrical 

results  as 

^  =  T^i)- J,.\  J^.\  ^  c  (3.6) 

Finally,  we  also  observe  that  midpoint  quadrature  rules  result  in 

I  u  -  Jx  ^  u(0  •  Jx( O 

J  u  ■  r/S  s:  u(C)  •  c>S(C')  (3.7) 

j*  u  u|C’)  • ' 

where  C  is  a  representative  center  point  of  the  figure. 

These  facts  summarize  the  principal  properties  of  normal  volume 
elements  that  we  shall  require. 

3.3.  The  Operators  div*  and  grad;^ 

The  fundamental  geometrical  properties  of  the  differential  operators 
div,  curl,  and  grad  stem  from  the  following;  div  relates  volume  integrals 
with  surface  integrals,  curl  relates  surface  integrals  with  line  integrals,  and 
grad  relates  line  integrals  with  endpoints.  We  will  now  describe  how  these 
properties  can  be  systematically  developed  for  the  div  and  grad  operators. 
div,,:  Gauss’  formula 

j  div  11  </I '  =  I  u  n  (IS 

when  applied  on  dV  using  the  quadrature  approximations  indicated  abo\e 
suggests  the  definition 

( 3.S ) 


r  div,,  u  —  ^  /1,u  bS 
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in  which  the  difference  operator  applies  to  opposite  face  center  points  . 
Using  the  component  representations  for  (he  terms  on  the  right, 

div,,  u  =  T  '  ^ /f,//' -  rXV,  (3.9) 

/ 

or,  noting  (3.6), 

div;,  u  =  (yco)  '  X  V  '''//1,-v  (3.10) 

/ 

grad,;.  The  approximations  to  be  identified  with  u  =  grad  ({>  will 
follow  from  the  line-integral  evaluation 

u  ■  d\  =  (f){B)  ~  <f){A)  (3.11) 

in  which  A,  B  are  end  points  of  a  curve  C.  We  will  describe  how  the 
covariant  components  of  u  with  respect  to  the  bases  [e(0)]  can  be  deter¬ 
mined  at  face  center  points. 

Tangential  Components.  Recalling  earlier  notations,  suppose  ). 

^y{Q\  +  )  ars  tangent  vectors  to  the  face  5S(0,  +  )  at  +  .  Evaluating  (3.1 1 ) 
along  the  line  connecting  the  centers  of  opposite  edges  /?i  ^.2+  and  /?,  ^  3  + 
we  obtain 

I  u-r/x  =  C(^i+.,+  )-C(^i  +  ,,  K  ,/=3.3 

-  C' 

in  which  we  have  identified  the  potential  variable  on  edges  of  dS  by  C. 
Next,  use  the  quadrature  evaluation 

u  ■  r/x  =  A,x  •  [u  •e,((2,  +  )]  =  A^x  ■  ^  ] 

•‘c 

in  which  u,  is  the  covariant  component  of  u(0|  *  )■  Thus, 

J  =  .M  ./  >'  7  =  -^  3 

=  A,i^{Qi,)  (3.12) 

Nontangential  Components.  We  will  also  require  evaluations  of  the 
non-tangential  components  of  u  at  the  center  points  of  faces  and,  as  for  the 
one-dimensional  case,  want  to  evaluate  these  in  terms  of  the  value  of  1/  at 
the  center  F  of  ST  and  of  (p  at  the  center  points  Q  of  its  faces.  Let  us  first 
introduce  the  one-sided  difference  operators 


a:(P(q,.)=  i-'/zi/^)] 


13.13) 
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Then  (3.8)  leads  to 

[  u- ^(0,  , ) 

in  which  the  line  integrals  are  to  be  taken  between  Q,  and  P  and  between 
P  and  respectively.  We  then  evaluate  the  integrals  by  one-sided 

quadrature  approximations  using  values  of  u  at  the  faces  to  obtain 

j  u  c/x  ^  (zJ,.v(P)/2)  )-e,((2,,  ) 

=  (A,x{P)ll)u,{Q,,) 

so  that 

[zi,.v(/>)/2]  ■  j )  =  ±vm.. )  -  4'(P)']  (3.14) 

We  can  summarize  this  construction  as  follows:  grad,,:  u(Q)  =  gvdd,,  (p 
expresses  the  set  of  relationships  between  the  covariant  components  of 
u{Q)  and  values  of  p,  i//,  and  C  at  points  Q  on  the  faces  of  an  element  given 
by: 

Tangential  Components: 

A,x-u^[Q)  =  A/e{Q)  (3.12) 

Nontangential  Components: 

{A,xl2)u,{Q)  =  Ap^m)  (3.14) 

where 

A-m)^  ±vm)-HP)^  (3.13) 

Equations  (3.14)  should  be  compared  to  Eqs.  (2.6).  By  adding  and 
subtracting  these  terms  we  obtain  an  equivalent  form: 

A,^{P)  P,u,[P)  =  A,(t){P) 

(3.15) 

iij{P)^H,(l>(P)~A,x{P)A,u,(P),A 


in  which  only  centered  operators  at  P  arc  involved  and  which  may  be  seen 
to  be  second-order  accurate  with  respect  to  P.  These  forms  play  an  impor¬ 
tant  role  in  establishing  energy  results. 


3.4.  A  Compact  Scheme 

We  will  now  return  to  the  primary  problem  of  solving  the  boundary 
value  problem  when  V  is  a  general  volume  w  ith  boundary  5.  Follow  ing  our 
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earlier  discussion,  we  describe  a  compact  scheme  for  solving  the  problem 
as  follows: 

(i)  In  each  element  solve 

div/,u=/,  u  =  grad;>  (3.16) 

Note  that  div/,  u  =  /  provides  one  equation  for  the  six  contravariant  com¬ 
ponents  of  u  on  the  faces  of  the  element,  while  u  =  grad/,  ^  expresses  three 
equations  for  the  covariant  components  of  u  on  each  of  the  six  sides  of 
in  which  the  tangential  components  are  evaluated  in  terms  of  the  variable 
C  at  the  center  points  of  edges  while  the  nontangential  components  are 
determined  by  the  variables  ^  and  ip. 

(ii)  Use 


=  (3.17) 

from  (3.4)  to  relate  the  contravariant  components  to  the  covariant  com¬ 
ponents  arising  in  (3.16). 

An  important  simplification  will  occur  in  this  construction  whenever 
e‘'  =  0,  iV  j.  When  this  situation  arises  we  call  the  scheme  strongly  com¬ 
pact;  otherwise  we  use  the  term  weakly  compact.  For  a  strongly  compact 
scheme  the  tangential  components  of  u  on  faces  will  not  affect  the  com¬ 
putation  directly  (and,  thus,  neither  will  the  variable  ()  but  may,  if 
required,  be  determined  in  terms  of  the  variable  C  by  a  postprocessing 
technique  using  the  completeness  relationships  described  below.  This  arises 
when  Cartesian  or  orthogonal  curvilinear  grids  are  involved. 

(iii)  Impose  continuity  conditions  for  u,  <p,  C  on  faces  common  to 
neighboring  elements.  In  case  is  a  boundary  element,  the  data  on  S  will 
be  assumed  to  determine  the  data  <p  and  II  on  any  face  tangent  to  S  by  the 
use  of  a  Taylor  series  approximation. 

We  are  required  to  calculate  u  and  <P  at  the  centers  of  the  faces  of 
elements,  ip  at  the  centers  of  the  elements  themselves,  and  (  on  the  edges 
of  the  faces.  By  examining  a  cube,  we  can  verify  the  fact  that  the  number 
of  unknowns  involved  in  the  system  of  equations  resulting  from  the  steps 
just  described  will  exceed  the  number  of  equations  by  just  the  number  of 
interior  edges  of  elements,  i.e.,  those  edges  that  are  incident  to  four  adja¬ 
cent  elements.  This  same  result  will  also  hold  for  a  more  general  covering 
{SV}.  In  order  to  make  the  scheme  determinate  we  will  add  the  following 
completeness  conditions. 

(iv)  Completeness  Conditions:  Let  N{R)  denote  the  center  points  P 
of  those  volume  elements  that  are  incident  to  an  edge  having  R  as  center 
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point.  Using  a  bilinear  interpolation,  express  ({R)  in  terms  of  the  values 
tl/{R)  which  are  associated  with  the  neighboring  elements.  The  result  can  be 
written 


(3.18) 

where  the  summation  includes  points  P  \n  N(R)  and  the  weights  w  are 
non-negative  and  sum  to  I.  As  noted  above,  when  the  scheme  is  strongly 
compact  (3.18)  will  allow  the  tangential  components  of  u  on  faces  to  be 
calculated  separately. 

3.5,  Energy 

We  will  now  show  how  an  energy  expression  can  be  derived.  From  it 
the  existence  and  uniqueness  of  the  solution  of  the  discrete  problem  will 
follow. 

Following  the  argument  given  for  the  one-dimensional  case,  use  the 
summation  rules  (2.9)  and  also  (3.18)  to  obtain 

T  div/,  ^  Ai(f)u‘  SSj 

i 

/ 

=  Z  [(•/'  +  5^/-^  <55,  +  A^xn,u,{n,u'  (55,)] 

I 

=  i/z  Z  A,u'  dS,  -f-  ^  A,.\fi,(u,u'  SS,) 

i  I 

or,  recalling  that  dS,  =  ^Je  A iX  A^x, 

T  div;,  =  (l)T  div/,  u  +  To  Z  J<-’) 

so  that  summing  over  volume  elements  yields 

ZdiV;,  (^UT  =  X</i<Jiv,,UT  +  Z  |z y(')|  To 

r  I'  r  I  (  3 

Finally,  using  the  telescoping  property  of  the  Zi  term  on  the  left-hand 
side  and  recalling  that  t  =  r,,  >/<'„,  we  find 

Z  ■  <5S  =  ^  ^  div;,  UT  -t-  ^  u  •  UT 

v  r  i' 


(3.19) 
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which  is  clearly  consistent  with  the  integral  form  of  the  energy  equation 
given  by  Eq.  (1.4).  Noting  that  (3.19)  will  ensure  that  (  is  bounded  in  norm 
by  ij/,  the  same  arguments  as  were  given  for  the  one-dimensional  problem 
will  imply  that  this  more  general  compact  scheme  converges  and  yields 
second-order  accuracy  for  i//,  (  and  u. 


3.6.  The  Potential  Formulation 

We  shall  indicate  how  the  scheme  just  described  can  be  formulated 
solely  in  terms  of  <}>,  tj/,  and  corresponding  to  the  potential  form  =  f. 
Again  following  our  discussion  of  the  one-dimensional  case,  we  write,  usinc 
(3.9)  and  (3.4), 


Tdiv/,  u  = 

I 

/  v  /  / 

or,  regrouping  terms: 

T  div^  u  =  X  +  Z  Z  ( -'  -0 ) 

/  /  \  ;  ^  /  / 

Using  (3.14)  in  the  first  summation  we  find 

Z  ^^e“u,  3S,  =  2  Z  A.e"  A ^cl>3S,l A, X 
/  / 

=  4  Z  3S,IA,x  -  mP)  Z  p.,c"  dSJA^x  ( 3.2 1  ) 

/  I 


in  which  the  summands  are  evaluated  at  the  centers  of  faces  about  the 
center  point  P  of  the  element.  The  second  group  of  terms  in  (3.20)  involve 
only  tangential  components  of  u  at  face  centers,  and  we  may  use  (3.12)  to 
express  these  in  terms  of  values  of  C  at  the  edges.  In  this  manner,  all  of  the 
terms  in  (3.20)  may  be  expressed  in  terms  of  the  values  of  <p,  i//,  and  2  in 
each  clement.  This  leads  to  the  potential  operator  div^,  grad,,  ({>. 

Since  «^,  ip,  and  C  arc  assumed  continuous  across  faces,  the  continuitv 
conditions  for  u  across  a  face  SS,  common  to  neighboring  elements  with 
centers  at  P  and  P'  reduce  to 

[w'(5A’,]  =0 

The  tangential  covariant  components  of  u  are  continuous,  since  they  are 


determined  by  Recalling  (3.14),  the  condition  that  the  nontangential 
covariant  component  m,  be  continuous  is  then 


(AfX  +  A,x')(l>  =  AjXtj/'  +  A,x'ip  (3.22) 

which  may  be  compared  to  (2.9a). 

Consider  the  treatment  ofV'<^  =  0  using  a  strongly  compact  scheme.  In 
this  case  =  and  the  condition  div,,  grad/,  =  0  as  expressed  by 

using  (3.20)  reduces  to 

^  H, €'•<!>  SS,/A,x  =  Pi,e“  dS,IA,x  (3.23 ) 

i  ^ 

which  may  be  written 

ij/ =  Y^<x,<f>  (3.24) 

I 

where 


(3.25) 


(3.26) 

we  can  conclude  that  the  maximum-  minimum  property  holds.  This  same 
conclusion  can  be  reached  for  a  weakly  compact  scheme. 

Referring  to  the  energy  expression,  it  will  be  seen  that  the  potential 
form  can  also  be  viewed  as  arising  from  a  variational  problem.  By  using  the 
completeness  relations,  the  computational  problem  will  thus  reduce  to 
determining  the  variables  (l>  and  tj/  associated  with  a  symmetric,  positive 
definite  matrix  and  thereby  allows  a  variety  of  familiar  solution  techniques 
to  be  used. 

3.7.  Degeneracy 

Recall  that  the  definition  of  a  normal  volume  element  allowed  for  the 
area  of  one  or  more  of  its  faces  to  vanish.  Such  degenerate  elements  play 
a  natural  role  at  the  boundary  5  and  are  especially  relevant,  as  will  be 
shown  shortly,  in  treating  problems  in  curvilinear  coordinate  at  the  origin. 


a,  =  (//,t’"  SSi/A,x)j(^  n^e"  3S,/A.x'^ 

I  =  1 

t 

If  we  write  the  continuity  condition  (3.22)  as 
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Let  US  suppose  that  a  face  with  center  Q  is  degenerate,  i.e.,  ^/S(0)  =  O. 
Referring  to  (3.8)  we  see  that  the  term  involving  the  factor  r/S(0)  will  not 
contribute  to  the  evaluation  of  div/,  u  and,  as  a  result,  the  contravariant 
component  of  u(0)  need  not  be  computed  on  the  face.  In  the  potential 
form,  this  simply  modifies  the  weight  a  associated  with  the  degenerate  face 
in  (3.25). 


3.8.  Curvilinear  Coordinates 


We  will  illustrate  how  the  discussion  applies  when  curvilinear  coor¬ 
dinates  are  used.  With  their  use  the  basis  (e^,  e^,  e^),  which  was  required  at 
each  face  of  a  volume  element,  can  be  constructed  analytically,  thereby 
considerably  reducing  the  preliminary  computational  steps  required.  Using 
the  potential  form,  we  shall  also  illustrate  the  treatment  of  the  situation  at 
the  origin  that  gives  rise  to  degeneracy. 

We  indicate  by  x  =  x(q)  a  mapping  in  general  curvilinear  coordinates. 
The  vectors 


dx 

dq' 


(3.27) 


are  tangent  vectors  on  the  coordinate  lines  formed  by  q'  =  const  and 
q'^  =  const  so  that  an  element  of  length  is  determined  by 

=  Z  Z 


in  which  are  metric  coefficients  given  by  g,y  =  g, -gy  Considering  [g]  = 
(g,,  gy,  g^)  as  a  covariant  basis,  a  reciprocal  basis  (g',  g\  g*)  is  given  as 

in  which  g  =  det(g,,).  Also,  if  g"  =  g'  gL  then  det(g''')  =  g  '■ 

Consider  a  curvilinear  volume  element  dV  which  is  formed  as  the 
image  of  a  rectilinear  volume  element  whose  volume  is  dq'  dq'  dq''.  Then 
the  edges  of  a  face  on  dV  given  as  ^/'  =  const  arc  coordinate  lines  on  which 
g^,  g;  are  tangent  vectors  and  the  oriented  area  of  the  face  is 

=  (g,  X  g  J  dq'  dq^ 

We  may  associate  a  normal  volume  element  51'  with  dV  by  construct¬ 
ing  tangent  planes  through  the  centers  of  the  surface  elements  through  dV. 
From  geometrical  considerations  it  is  evident  that  if  we  make  the  following 
identifications 


),.\c,  =  ^,dq\ 


/=  1.2,3 


(3.28) 
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the  results  of  our  earlier  discussion  will  apply.  If  i/'(g)  indicates  a  con- 
travariant  component  of  u  in  the  curvilinear  system  we  find,  in  place  of 
(3.10), 

divy,  u  =  '  Z  Jg  »'(g)A/^/'  (3.29) 

I 

Similarly,  u  =  grad,,  expressed  earlier  by  (3.12)  and  (3.14),  now  result  in 
the  following  components. 

Tangential  Components: 

^r/'-n,(g)  =  zl/  (3.30) 

Non  tangen  t  ial  C omponen  ts : 

(Jq'/2)u,(g)  =  ^r<l>  (3-31) 

Examples 

The  cases  of  cylindrical  and  spherical  coordinates  may  help  illustrate 
these  results. 

Cylindrical  Coordinates.  Here  q  =  (r,  0,  r)  and  a  calculation  gives 
Jg=^r  and  Jg  =  r,  Jg  g-  =  l/r,  Jg  g^^  =  r.  The  center  points  of  faces 
Qi+  have  the  coordinates 

Px±=ir±  ArjX  0,  z\  ?2±  =  (r,  ^  ±  ^0/2,  z\  Py,  =  (r,  0,  ±  Az/2) 

Adopting  more  standard  finite  difference  notations,  let  <j>{i  Ar,  j  AO,  k  Az)  = 
^rid  extend  the  meaning  of  the  spatial  averaging  and  differencing 
operators  p,  A  to  this  case  by  writing  (//,.,  pp,  p.),  (A^,  Ap,  A.).  We  find  that 
in  a  cell  with  center  at  T  =  (/  Ar,  j  AO,  k  Az) 

4  (  1 

div„  grad;,  =  -  {yJP  ^  i.z.j.k  +  ('’  -  ^r/2)]  0,^  ,  2.,.a- 

+  Oo4>,.,.k  +  (3  32) 

in  which 

Assuming  equal  spacings,  the  continuity  conditions  across  faces  have  the 
form 

4>  =  l4  (3.33) 


in  which  p  =  (p,,  pp,  pA. 
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A  cell  with  center  at  (Aril,  j  AO,  k  Az)  has  the  origin  as  a  degenerate 
face.  In  this  case  we  see  that  the  coefficient  of  the  value  of  the  unknown 
^0  ,  ^  vanishes,  so  that  no  assumption  other  than  boundedness  at  the  origin 
is  required  for  the  scheme. 

For  cells  not  involving  boundary  points,  we  may  use  the  continuity 
condition  in  (3.29)  to  express  the  result  solely  in  terms  of  tp.  By  adopting 
the  notation  S^  =  A  JAr,  etc.  to  indicate  divided  differences,  the  potential 
operator  for  such  cells  reduces  to  the  form 

div/,  grad/,  .^=:  !<:>;  + +  ^  + ip  (3.34) 

which  can  be  seen  to  be  consistent  with  the  familiar  differential  form  in 
cylindrical  coordinates.  However,  we  emphasize  again  that  the  boundary 
conditions  involving  (p  cannot  be  handled  directly  from  this  form  but  must, 
instead,  be  developed  using  (3.30)  directly. 

Spherical  Coordinates.  In  this  case  q  =  {r,(p,0)  and  a  calculation 
shows  g  =  r^  sin  (p  while  y/gg''=r^,  ^/g  g^^  =  sin  tp,  >Jgg^^=\.  Again 
assuming  a  uniform  mesh  spacing,  we  find  in 

4  r  sin  y 

div/,  grad/,  ^  +  laj.u  +  {r-Arjiy- 

+  ^'•>+  '/2.A-  +  sin(z  -  AyJ2)(p,_  ,  ^  J 

X  -  Oo'A,.y.*|  (3.35) 

in  which 


(T„  =  {r~  -f  Ar^ll)  sin  x  +  -^  (sin  i  cos  AyJ2)  -f 

[Ar  Ax  A0~) 

With  //  =  (/<,., /i^, ////)  the  continuity  conditions  again  lead  to  <p  =  p\p. 
We  see  that  the  situation  at  the  origin  is  exactly  similar  to  that  in  the 
previous  example;  in  particular,  only  the  boundedness  of  the  solution  at 
the  origin  is  required.  Away  from  boundary  cells,  the  continuity  conditions 
can  be  used  to  reduce  matters  to  the  form 

-4,1^  +  Ar^  12)  sin  <5^  +  2  sin  / 

+  {s\nycosAxl2)6]  +  dl\p  (3.36) 
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which  can  be  seen  to  be  consistent  with  the  potential  form  of  the  diffusion 
equation  in  spherical  coordinates.  We  again  add  a  precautionary  note 
about  attempting  to  treat  the  boundary  conditions  involving  <f>  from  this 
form. 

3.9.  The  curl^  Operator 

In  our  treatment  of  u  =  grad,,  <f),  emphasis  was  placed  on  evaluating 
the  tangential  components  of  u  at  the  center  Q  of  each  face  (5S  using  a 
covariant  basis  [ef^)].  Using  (3.12),  these  were  e.xpressed  in  terms  of  the 
variable  c  which  was  associated  with  the  center  points  of  the  edges  of  the 
face.  The  completeness  conditions  (3.18)  were  then  used  to  relate  (  to  i//. 
The  contravariant  component  of  u  which  arose  in  the  definition  of  div,,  u 
were  evaluated  in  terms  of  these  covariant  components  using  (3.4). 

It  is  desirable,  also,  to  also  have  a  finite  volume  approximation  to  the 
curl  operator,  which  we  will  indicate  by  curl,,  u.  A  satisfactory  calculus 
would  then  yield  both  div,,  curl,,  u  =  0  as  well  as  curl,,  grad,,  =  0  as  iden¬ 
tities.  We  shall  now  indicate  how  this  may  be  accomplished. 

In  Fig.  2,  T  indicated  a  representative  vertex  of  an  element.  Let  [e(r)] 
indicate  a  basis  formed  by  the  edges  which  intersect  at  T  using  the  orienta¬ 
tion  established  by  [e(P)];  also  let  e{R)  indicate  the  translation  of  e{T)  to 
the  center  point  R  of  an  edge  TT'.  A  suitable  definition  of  curl,,  is  suggested 
by  the  following;  If  cu  =  curl  u,  then  a  quadrature  approximation  in  Stokes’ 
formula 

r  r 

(D(IS=  ur/s 

leads  to 

co  SS^A^.x-  Aj(\\  •  )  -  zl^.v  •  •  e,) 

involving  the  edges  of  a  face  ^S.  Using  component  representations,  we  are 
led  to  define  curl,,  in  terms  of  the  set  of  contravariant  components  to' 
associated  with  the  faces  of  dV  hy  means  of 

0}  =  curl,,  u  ou'  =  (AfiJA^x  -  liJA^x)  (3.37) 

in  which  the  covariant  components  li^u  c.  With  this  definition  the 
identity 

div,,  curl,,  u  =  0  (3.38) 

follows  by  noting  the  cancellation  of  “line  integrals”  along  the  common 
edge  of  faces.  [It  is  also  possible  to  identify  a  vector  associated  with  the 
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center  point  of  dV,  say  (curl/,  u>  =  X! but  this  vector  will  not  be 
annihilated  by  div,,  and  for  that  reason  will  not  be  considered  further 
here.  ] 

As  already  noted,  our  earlier  treatment  of  u  =  grad,,  <f)  will  not  furnish 
a  direct  means  of  evaluating  the  covariant  components  u  occurring  above 
and  thus  we  cannot  interpret  curl/,  grad/,  In  order  to  overcome  this, 
introduce  a  new  variable  x  which  will  be  defined  at  vertices  T  and  is  often 
called  a  box-variable.  Recalling  the  intrinsic  definition  of  grad  indicated  by 
(3.1 1 ),  we  will  define  grad/,  along  edges  by 

A,x-u,(R)  =  Arx(R)  (3.39) 

in  which  AjX  is  the  distance  between  vertices  along  an  edge  with  direction 
e,  on  which  R  is  the  center  point  and  A,x{R)  is  the  difference  in  values  of 
X  at  the  vertex  neighbors  of  R.  Then 

co'(u)  =  {AjAkX  -  ^  =  0 


i.e., 


curl/,  grad/,  =  0  (3.40) 

This  suggests  a  theoretical  advantage  in  using  the  box-variable  x. 
instead  of  C  which  was  used  in  the  earlier  discussions.  The  most  direct  way 
to  accomplish  this  is  to  define  the  edge  variable  ((R)  as  the  average  of  /. 
at  vertices  neighboring  R,  i.e.,  (  =  and  also  reexpress  the  consistency 
conditions  (3.18)  so  that  now  xXT),  rather  than  C(^)»  will  be  related  to 
il/iP)  at  each  interior  vertex  by  a  bilinear  interpolation  of  the  form 
X  =  Qiil/),  which  will  now  involve  the  six  neighboring  elements  having  T  as 
a  vertex. 

4.  THE  TIME-DEPENDENT  PROBLEM 
4.1.  The  Potential  Form 

The  potential  form  of  the  time-dependent  diffusion  equation  (1.1)  is, 
with  /  =  0, 


4> ,  =  i\\\  gxad  (!>  (4.1) 

The  discrete  finite  volume  operator  div/, grad/,0  corresponding  to  divgrad0 
has  been  shown  to  yield  a  discrete  energy  expression  which  corresponds 
to  the  time-independent  energy  terms  in  (3.19).  The  principal  question 
remaining,  therefore,  is  how  to  include  the  term  <f>,  in  the  numerical  scheme 
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which  will  provide  the  approximation  to  the  term  {dldi)\  <f>^  dV  in  the 
energy  expression. 

Divide  [0,  T]  into  N  equal  parts  of  width  At,  write  t„  =  nAt,  and 
adopt  the  standard  finite  difference  notation  u{n  At)  =  u".  Also,  let  fi,,A, 
indicate  the  central  averaging  and  difference  operators  which  effect  the  time 
index  as  indieated; 


From  (2.4)  it  is  evident  that 

We  will  consider  the  following  three-level  compact  scheme. 
A.  Compact  Scheme. 


(4.2) 


(4.3) 


A/ij/'' =  At  d\v,,gT2L(i,,(f>"  (4.4) 

=  (4.5) 

to  which  the  appropriate  continuity  and  completeness  conditions  (3.22) 
and  (3.18)  are  understood  to  apply.  By  multiplying  (4.4)  by  ij/"  we  obtain 

=  At  ■  ij/"  div,,  grad;,  <!>"  (4.6) 

in  a  volume  element.  The  term  on  the  right-hand  side  has  already  been 
shown  to  lead  to  the  energy  expression  (3.19).  Thus,  by  summing  (4.6)  over 
elements  in  a  time-strip  and  using  the  implied  continuity  conditions  across 
faces  we  will  obtain 

X  1  +  AtY^n"  ■W'x  =  AtY^  ■  dS  (4.7) 

y  y  s 

which  is  the  discrete  energy  expression  that  was  expected. 

Assuming  dissipative  boundary  conditions  in  the  sense  of  Kreiss 
(1964),  it  is  evident  from  this  result  that  we  can  conclude:  both  ip"  and  (p” 
are  hounded  in  a  discrete  norm  by  the  initial  and  boundary  data.  Because  the 
scheme  is  consistent,  this  also  implies  that  the  scheme  converges.  (Our 
previous  discussion  of  the  maximum  minimum  property  can  allow  us  to 
conclude  convergence  as  well.) 

We  can  summarize  this  in  several  equivalent  ways,  in  each  of  which 
the  additional  continuity  and  completeness  conditions  are  assumed  to 
apply; 
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B.  Two-Step. 


and 

+  1  _  +  1/2  ^  div,,  grad,  (4"  +  (4  8b) 

As  above,  this  is  a  two-step  scheme  in  which  the  first  requires  solving  the 
equations  in  a  time  strip  for  and  implicitly  in  terms  of  the 

boundary  data  for  ^^(/„  +  i/2)  and  the  initial  data  ip".  With  these  values 
determined,  then  ip"^'  can  be  determined  explicitly  from  the  second  set  of 
equations  and  furnishes  initial  data  for  the  next  time  strip. 

C.  One-Step. 

Eliminating  ip"  in  the  above  scheme  we  find 

=  J/^,div,  grad,  (4.9) 

This  is  a  one-step  Crank-Nicholson  form. 

Example.  Using  our  earlier  one-dimensional  notations,  consider  the 
treatment  of  (l),  =  (p^^.  We  find  that  the  one-step  scheme  may  be  written  in 
potential  form  as 

A ,\p'- =  Kp,{p(p'l —  \p1),  /e  4  (center  points)  (4.10a) 

where  K  —  Atjh^,  h  =  Axl2.  To  this  we  add  the  continuity  conditions  in  the 
form 


— =  /e  4  (end  points)  (4.10b) 

Comparing  with  (2.9)  we  again  recognize  that  the  variables  (p,,  (/',  are  the 
odd-even  components  of  a  tridiagonal  system  and  are  therefore  readily 
obtained  at  each  time  step. 

Numerical  tests  using  the  solution  ^  =  exp(/  — .v)  with  initial  condi¬ 
tions  at  /  =  0  produced  the  following  (max)  error  norms  for  (p  and  (p^  at 
/  =  1  in  which  At  =  Ax: 


Total  number  of  points 

(p  error 

(p^  error 

9 

0.00116002 

0.0186447 

17 

0.00029997 

0.00485956 

33 

0.0000768122 

0.00123751 

These  results  appear  to  verify  the  accuracy  asserted  by  the  theory. 
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4.2.  An  ADI  Scheme 


Except  for  the  one-dimensional  example  just  considered,  the  basic 
scheme  (4.4)  is  implicit  and  effective  solution  methods  must  be  considered 
in  order  to  treat  the  general  problem  in  higher  dimensions  in  a  practical 
manner.  It  is  especially  desirable  that  a  proposed  method  also  provide  an 
effective  means  of  treating  the  steady-state  problem  as  well.  We  will  now 
show  that  a  Peaceman-Rachford-type  ADI  scheme,  using  a  straightforward 
treatment  of  intermediate  boundary  conditions,  can  solve  the  finite-volume 
problem  with  second-order  accuracy  in  both  space  and  time.  Again,  the 
existence  of  an  energy  estimate  will  provide  the  key  to  convergence. 

We  will  find  it  convenient  to  make  a  few  slight  changes  in  some  of  our 
notations.  First,  we  will  let  v  indicate  the  time  average  of  a  variable  v,  i.e., 
let 

=  (4.11) 


Also,  we  let 


F"  =  t-'J,(u''-^S,) 


(4.12) 


so  that  (3.9)  leads  to 

div,,  u”  =  F" -h  -h  F; 

Then,  for  example,  the  one-step  scheme  (4.10)  can  be  written  as 

d,i//''  =  dt(F';  +  F'^-hF';)  (4.13) 

in  which  u  =  grad„  (f>.  Finally,  we  indicate  by  g"  assumed  boundary  data  for 
on  5"  at  time 

Limiting  our  discussion  to  two  dimensions,  consider  the  following 
two-step  scheme; 


=  {At/2){F'\ "  +  F")  (4.14a) 

>  _^"+'/2  =  (j^/2)(F';"  ’/-VFr ’)  (4.14b) 

and  note  its  similarity  to  (4.8).  Recalling  our  discussion  of  the  one-dimen¬ 
sional  example,  it  is  evident  that  the  first  step  can  be  solved  as  a  one- 
dimensional  problem  using  boundary  conditions  -  =  g"^  ’  ■.  We  solve 
the  second  step  similarly,  but  with  the  boundary  conditions  <f)''  *  '  —  g"  *  '. 

Adding  and  subtracting  (4.14a)  and  (4.14b),  and  relabeling  the  time 
level  for  purposes  of  later  discussion,  we  find 

=  +  F'^) 

if/"  -ij/"  =  {Al/4)(F]  '^^-FV  ''') 


(4.14a) 

(4.15b) 
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The  consistency  of  (a)  with  the  one-step  Crank-Nicholson  form  (4.11)  is 
evident,  the  truncation  error  being  O(At^).  Equation  (4.15b)  is  also  consis¬ 
tent  with  (4.5)  since  the  right-hand  term  is  O(At^).  We  may  regard  (4.14) 
either  as  a  scheme  to  approximate  the  solution  of  the  compact  scheme  (4.4) 
or  as  an  independent  scheme  to  solve  the  differential  equation. 

We  shall  show  that  this  scheme  leads  to  an  energy  estimate  which  is 
consistent  with  (4.7)  to  within  terms  of  0{At^).  Recalling  the  basic  energy 
argument  in  (2.11)  and  the  definition  of  F,  in  (4.12)  we  find 

ij/  ■  F^  —  A{<l>u'  dS^ )  —  ^|(iFm,  ) 

yl)  ■  F2  =  A{^u^  dSj)- 

so  that  by  multiplying  (4.15a)  by  ^  and  then  using  (4.15b)  -  e  find 

:^A,{i(/'’)^  X  +  [/i,(M'w, )"  -I-  a^2(w^W2)"] 

2  At 

=  lA{(j)u^  dS^ )  +  A{^rd  r/52)]  +  0{At") 

Summing  over  volume  elements  and  noting  that  if/  =  if/  +  OiAr),  etc.,  we 
again  obtain  the  energy  expression  (4.7)  to  within  terms  O(At^).  This  result 
is  sufficient  to  enable  us  to  conclude  that  the  ADI  scheme  converges  and 
furnishes  an  approximation  with  the  same  degree  of  accuracy  as  the  com¬ 
pact  scheme  (4.4).  This  result  is  independent  of  any  of  the  fixed  ratios 
At/{Ax,y  occurring  in  either  scheme. 

To  treat  the  three-dimensional  problem,  consider  in  place  of  (3.14a), 
etc.,  the  following: 

^|^n  +  1/4  _  ^  {At/4){2F’;  +  F"  +  '/“’) 

+  1/2  _  ^  (^//4)f /T"  4  1/4  +  2F"  ^  '/') 

'  ^  (4.16) 

^  3/4  _  +  1/2  ^  ^  -3/4  ^  2F"  ^  '/^) 

+  1  _  +  .1/4  _  (^;/4)(2F';  ^  +  F"  ^  ^/‘’) 


to  which  the  intermediate  boundary  conditions  =  g"  ^  k  =  1/4,  1/2, 
3/4,  I  apply.  Then 


^  =  AtliF'l  ^ '  -h  F")/2  -t-  (F''  * 


.1'4 


+  F"  ^ 


1/4 


)/2  +  F 


t  1/2  ■ 


so  that  the  .scheme  is  consistent  with  the  Crank-Nicholson  form.  Also, 

(lA"  " '  +  ij/”)/!  -  lA"  *  =  (/1//4)[(F';+  '  -  F")  +  (F"2  *  -  F’l  ^  '  ")] 

The  last  pair  of  equations,  which  correspond  to  (4.15)  in  the  two-dimen¬ 
sional  case,  can  provide  the  basis  for  an  energy  argument,  although  we 
omit  the  details  here. 
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5.  CONCLUSIONS 

We  have  described  a  finite  volume  method  which  closely  maintains  the 
parallel  between  differential  and  diflerence  arguments.  By  using  intrinsic 
geometrical  properties  of  the  elements,  we  were  able  to  describe  versions  of 
the  div,  curl,  and  grad  operators,  which  led,  using  formal  summation  by 
parts  techniques,  to  discrete  energy  expressions  as  well  as  to  the  identities 
div/,  curl,,  u  =  0  and  curl,,  grad,,  ^  =  0.  The  solution  of  the  initial-boundary- 
value  problem  for  the  diffusion  equation  was  described  directly  in  terms  of 
these  operators  by  compact  schemes  and  the  resulting  energy  equations 
ensured  convergence.  The  schemes  could  also  be  simplified  to  a  potential 
form  which  can  offer  computational  advantages,  and  the  treatment  of 
general  curvilinear  coordinates  was  shown  to  result  from  a  specialization  of 
these  results.  Numerical  examples  validated  the  second-order  accuracy  of 
both  (f)  and  grad,,  <f)  which  the  theory  predicts. 
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An  Implicit  Enthalpy  Scheme  for  One-phase  Stefan 
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Milton  E.  Rose  t 


Introduction. 

When  finite  difference  equations  can  be  used,  ADI  schemes  provide  an 
effective  numerical  method  to  solve  multidimensional  diffusion  equations  because 
they  require  the  solution  of  simple  tridiagonal  systems  of  equations 
corresponding  to  a  one-dimensional  operator  at  each  time  step.  Although  such 
schemes  can  be  expected  to  have  similar  advantages  for  the  treatment  of  Stefan- 
type  problems,  few  studies  of  their  application  to  this  class  of  problems  have 
appeared  in  the  literature.  Several  studies  by  White  [5,6]  have  suggested  that  the 
treatment  of  the  underlying  nonlinear  one-dimensional  operator  which  arises  in 
such  applications  may  require  a  significant  number  of  iterations  to  converge  and 
may  also  require  a  limitation  in  the  time  step,  thus  offering  only  a  moderate 
improvement  over  explicit  or  SOR-type  iterative  methods.  This  difficulty 
apparently  ari  ;es  because  standard  finite  difference  operators  may  not  produce  a 
symmetric,  positive-definite  operator  for  the  Stefan  problem. 


*  Tliis  research  was  sponsored  by  the  Air  Force  Office  of  Scientific  Research 
under  Contract  No.  AFOSR  F49620-89-C-0010. 

1"  Present  address:  4505  Tower  Rd.  Greensboro,  NC  27410 
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This  paper  examines  a  novel  enthalpy  formulation  of  one-phase  Stefan 
problems  which  is  suggested  by  a  more  general  compact  finite-volume  treatment 
of  diffusion  equations  (Rose[4]).  Limiting  our  present  study  to  one-dimension, 
we  introduce  a  discrete  operator  which  is  symmetric  and,  we  conjecture, 
positive-definite  and  can  be  solved  on  a  fixed  solution  domain  with  a  relatively 
few  iterations.  The  numerical  results  indicate  that  the  position  of  the  melting 
front  can  be  determined  with  first  order  accuracy  by  this  method,  the  number  of 
iterations  at  each  time  step  being  determined  largely  by  the  number  of  cells 
traversed  by  die  front  during  a  time  step. 

We  first  review  the  enthalpy  method  and  a  standard  finite  difference 
scheme  usually  associated  with  the  diffusion  operator.  We  then  propose  a 
compact  finite  difference  scheme  for  solving  a  one-phase  Stefan  problem  and 
conclude  by  presenting  numerical  results  for  several  model  problems. 
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1  Enthalpy  method. 

The  one-dimensional  Stefan  problem  involving  the  enthalpy  e,  heat  flux 
q,  and  temperature  T  is  described  by  the  following  parabolic  system: 


(a) 

et  +  qx  =  0 

(conservation  of  energy) 

(1) 

(b) 

q  =  -KTx 

(Fourier's  law) 

(c) 

T  =  T(e) 

(equation  of  state). 

in  which  k  is  the  coefficient  of  tliermal  conductivity  and  a  constant  unit  density  is 
assumed.  For  ordinary  heat  conduction  problems  the  equation  of  state  is  a  one- 
to-one  transformation,  in  which  case  e  can  be  eliminated  and  the  familiar  heat 
conduction  equation  for  T  results.  When  melting  can  occur  (involving  a  latent 
heat  L)  the  equation  of  state  can  be  described  by  the  model  form 

T-Tq  =  e,  0  <  e  (2) 

=  0,  -L<e<0 

=  y-  (e  -  L),  e  <  -  L 

in  which  1/  are  specific  heats  and  Tq  is  a  a  reference  temperature  which  we 
may  assume  vanishes,  i.e.,  Tq  =0.  (We  refer  the  reader  to  Crank  [1]  for  a  fuller 
discussion). 

In  the  problem  to  be  discussed  below,  we  will  seek  to  determine  (e,q,T) 
when  a  slab  of  material  (such  as  ice)  has  the  initial  value  e  =  -L  and  is 
subsequently  melted  by  the  addition  of  heat  through  one  or  both  of  its  faces  (at 
the  endpoints,  say,  of  the  unit  interval  0  <  x  <  1)  as  specified  by  a  linear 
combination  of  the  temperature  and  flux.  Of  particular  interest  is  the  position  x 
=  s(t)  of  the  interface  between  the  solid  and  liquid  phases.  By  integrating  (1) 
across  this  interface  we  obtain  with  the  use  of  Gauss’  theorem 

(e]  s'(t)  =  [q]  ,  (3a) 

[T]  =0, 

in  which  [  ]  indicates  the  jump  in  value  across  x  =  s(t).  Across  the  interface,  [e]  = 
L  and  we  thus  obtain  the  familiar  free  boundary  conditions 

L  s'(t)  =  [q],  (3b) 

T  =  0  ,  X  =  s(t) 
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which  were  first  described  by  Stefan. 

Because  (3)  is  obtained  by  integrating  (1)  over  a  domain  including  the 
interface,  the  solution  (e,q,T)  can  be  shown  to  be  a  weak-solution  of  (1)  which  is 
smooth  in  each  subdomain  separated  by  the  interface.  Following  a  suggestion  of 
Lax  concerning  conservation  laws,  such  solutions  can  be  expected  to  result  as  the 
weak  limits  of  appropriate  discretized  forms  of  (1).  In  [3],  Rose  formulated  the 
Stefan  problem  in  the  manner  described  above  fenthalpy  methodl  and  examined 
the  use  of  a  simple  explicit  numerical  scheme  to  obtain  the  solution.  (A  fuller 
discussion  of  this  method  may  be  found  in  Crank  [1]). 

Let 


E  =  1/2  d/dt 


J 


1 

Tedx  + 


1 


1 

^Tqdx  +  Tq 


(4) 


represent  the  total  "energy"  on  the  interval.  By  multiplying  (la)  by  T, 
integrating  over  the  fixed  interval  [0,1]  taking  account  of  the  energy  balance 
condition  (3b)  across  an  interface,  and  noting  the  assumption  that  Tq  =0  we 

arrive  at  an  "energy  balance"  equation  with  the  following  result: 

E  =  To[q]=0.  (5) 


In  discussing  smooth  solutions  to  the  diffusion  equation.  Rose  [4]  consU'ucted  a 
finite  volume  scheme  which  was  consistent  with  (1)  and  also  satisfied  a  discrete 
version  of  (5).  As  a  result,  the  scheme  could  be  shown  to  converge  and  produce 
second  order  accuracy  in  both  T  and  q.  By  adapting  this  scheme  to  the  Stefan 
problem  we  expect  that  an  energy  result  similar  to  (5)  will  also  emerge,  thereby 
resulting  in  a  symmetric,  positive-definite  operator  which  is  particularly  relevant 
to  the  construction  of  a  weak  solution  to  the  problem.  (Because  this  scheme  can 
also  be  interpreted  as  a  non-conforming  finite  element  scheme,  the  error  analysis 
given  by  Elliott  [2]  for  the  enthalpy  method  may  not  be  entirely  applicable.) 
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2  A  standard  finite  difference  scheme. 

Consider  the  two-dimensional  diffusion  equation 

Ut  =  Uxx  +  Uyy  (4) 

and  let  F,  G  denote  standard  finite  difference  approximations  to  Uxx  and  Uyy 
respectively  both  involving  a  spatial  parameter  Ax.  Using  standard  finite 
difference  notations,  an  ADI  solution  technique  to  solve  (4)  is 

(a)  (u  u  ^)/T  +  (5) 

(b)  (u  -i-G 

where  t  =  At  /  2.  This  scheme,  which  is  formally  second  order  accurate  in  both 
At  and  the  implied  spatial  parameter  Ax,  is  unconditionally  stable.  The 
effectiveness  of  this  scheme  arises,  in  part,  from  the  observation  that  each  step 
involves  the  solution  of  a  tri-diagonal  matrix  whose  solution  can  be  readily 
obtained  (Thomas  algorithm).  In  one  dimension  this  scheme  reduces  to  the  form 


(a) 

(u 

_  p  n+1/2 

(b) 

(u 

_  p  n+1/2 

so  that  in  this  case  only  the  first  step  (a)  requires  the  solution  of  an  implicit 
system  of  equations. 

For  a  model  one-phase  problem,  e(x,0)  =  -L  initially  and  we  may 
assume  -  L  <  e(x,t)  for  t  ^  0.  If  we  introduce  the  notation 

y(e)  =  y'^  ,  0  <  e  (7a) 

=  0,  -  L  <  e  <0 

the  equation  of  state  (2)  can  be  written 
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T(e)=  Y(e)e,  -L<e  (7b) 

which  is  the  form  we  shall  assume  in  the  following  discussion. 

Corresponding  to  (6),  a  conventional  difference  scheme  arising  from 

(1)  is 

(a)  (ei  +  (qi+i/2  -  qi_i/2  /  Ax  =  0 

(8) 

(b)  (q  n+1  -  ei  "+ 1«)  / 1  +  (qj+i/i  “+>«  -  qi_i/2  "+'«)  /  Ax  =  0 

with 


qi. - Ki(Ti+i/2-Ti_i/2)/Ax.  (9) 

in  which  T  =  T(e)  is  given  by  (6).  As  a  result,  F  has  the  interpretation 

Fi  =  [Ki+i/2(Ti+i  -  Ti)  -  Ki_i/2(Ti  -  Ti_i)]  /  Ax  2  (10) 

=  A(KATi) 


which  is  the  standard  2nd  order  difference  operator  corresponding  to  (kux)x» 

Here  the  subscript  i  has  the  range  [1,M— 1],  corresponding  to  a  subdivision  of  the 
interval  (0,1)  into  M  equal  subintervals  each  of  length  Ax  =  1/M.  Thus,  for  n  = 
0,  1,  2  ...  (8a)  leads  to  a  nonlinear  system  of  tridiagonal  difference  equations 

ejn+l/2  _  X  A  (KAT(ei"+^/2  ) )  ^  ^.n 

With  so  determined,  ej  can  be  obtained  explicitly  from 

ei'^+^  =  T  A  (KAT(ei"+*/2  ) )  ^  (115) 


As  noted  earlier,  the  explicit  scheme  defined  by 
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ei  =  T  A  (KAT(ei“  ))  +  ei*'  ,  n  =  0,1,2,... 

was  originally  used  in  [3]  to  demonstrate  the  feasibility  of  the  enthalpy  method; 
although  it  requires  the  severe  stability  condition  At  /  Ax^  <  1/2  it  suggests  that 
the  implicit  scheme  described  by  (11)  may  form  the  basis  for  an  approach  to 
solving  the  Stefan  problem  without  this  stability  condition.  However,  we  were 
not  able  to  develop  a  satisfactory  method  to  solve  this  system  effectively  (c.f. 
White  [5,6]),  possibly  because  the  operator  in  (11a)  may  fail  to  be  symmetric, 
positive-definite  when  (7b)  holds  .  We  will  now  describe  a  slightly  modified 
descritization  method  which  appears  to  allow  for  an  effective  solution  technique. 


3  A  compact  finite  difference  scheme. 

If  we  examine  the  conventional  scheme  represented  by  (7),  (8),  and  (9) 
we  observe  that  (8)  is  consistent  with  the  conservation  equation  (la)  on  intervals 
{  xj  -  Ax/2  <  Xj  <  xj  +  Ax/2}  while  (9)  has  been  made  consistent  with  Fourier's 

law  (lb)  by  relating  the  flux  at  the  endpoints  of  each  such  interval  to  the 
temperature  values  Tj  Tj+|  related  to  the  centerpoints  of  an  interval  and  its 

neighbors.  In  particular,  notice  that  (9)  will  provide  second  order  accuracy  only 
if  the  coefficient  k  is  smooth  across  subintervals. 

We  will  now  describe  a  slightly  different  numerical  scheme  which  is 
suggested  by  a  domain-decomposition  argument.  With  a  fixed  choice  of  disjoint 
subintervals,  we  propose  to  describe  the  physics  within  each  subinterval  by 
solving  (1)  in  isolation  from  its  neighboring  intervals  (i.e..  compactly)  and  then, 
separately,  extend  the  solution  in  the  large  by  means  of  continuity  conditions  for 
e,  q,  and  T  across  the  boundaries  of  contiguous  subintervals.  For  sufficiently 
small  subintervals,  an  approximate  solution  of  (1)  in  each  subinterval  which 
satisfies  initial  and  boundary  conditions  can  be  obtained  by  use  of  a  Taylor  series 
With  the  addition  of  the  remaining  continuity  requirements,  the  result  will  lead, 
to  a  compact  discrete  approximation  of  the  problem  throughout  the  fundamental 
domain. 
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In  order  to  implement  this  idea  for  the  one-dimensional  problem, 
introduce  the  M  disjoint  intervals  Jj  =  {  Xj- Ax  <  xj  <  Xj  +  Ax},  i  =  1,  3,  ..., 

2M  -  1  as  suggested  by  the  stencil 

1  •  1  •  I  •  I 

Xj-2  Xf  Xi+2- 

in  which  the  centerpoints  have  been  indicated  by  • .  First  apply  the  conservation 
law  (la)  in  each  interval  I[  to  obtain  the  approximate  condition 

(a)  (ei  -  ei  /  T  -h  (ql+i  _  q+._^  n+1/2)  /  (2^^)  =  0 

(12) 

(b)  (ei  -  ei  j  ^  ^  n+1/2  _  q+._^  n-H/2)  /  (2Ax  )  =  0 

i=  1,3,...,  2M-1 

involving  the  right  and  left  limits  q-  of  the  fluxes  at  the  endpoints  of  the  interval 
Ji  (compare  to  (8) ).  Next,  temporarily  omitting  the  superscripts,  express 

Fourier's  law  (lb)  for  these  flux  values  by  backward  and  forward  differences  in 
the  following  manner: 

q-i+ 1  =  -  K-i+ 1  (T-i+ 1  -  T(ei) )/  Ax  ,  (13) 

qVl  =  -kVi  (T(ei)-TVl  )/Ax. 

in  which  T(ei)  is  given  by  the  equation  of  state  (7b),  viz., 

T(ei)  =  7(ei)  ej  ,  -  L  <  ei  (7b) 

Clearly,  (7b),  (12),  and  (13)  are  consistent  with  the  physical  system  (1)  in  each 
subinterval.  In  addition,  we  will  impose  the  following  continuity  conditions 
across  the  endpoints  of  intervals: 
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[q]i  =  [T]i  =  0 


i  =  2,  4, ....  2M-2  . 


(14) 


Because  each  interval  is  fixed  in  time,  these  conditions  are  consistent  with  tlie 
general  application  of  Gauss'  theorem  across  interval  endpoints  as  expressed  by 
(3a). 


Following  [3],  we  call  the  scheme  described  by  (7b),  (12),  (13).  and 
(14)  a  compact  finite  difference  scheme  to  solve  (1).  Note  that  the  variables 
which  occur  are:  (T^,  qi)  at  cell  endpoints  and  e  at  cell  centerpoints.  The 

number  of  variables  can  be  reduced  by  setting  q  =  q"*"  =  q~  and  also  T  =  =  T“, 

so  that  the  continuity  conditions  [q]  =  [T]  =  0  will  automatically  be  satisfied. 

Thus,  except  for  boundary  values,  the  primary  variables  can  be  reduced  to  the 
pair  (T,  q)  at  cell  endpoints  and  e  at  cell  centerpoints.  (For  the  case  of  simple 
heat  conduction  (y(e)  =  1)  it  is  a  simple  matter  to  show  that  the  system  of 
algebraic  equations  resulting  from  (7),  (12),  (13),  and  (14)  is  uniquely 
determined  by  the  initial  and  boundary  conditions  for  the  discrete  problem) 

Although  the  definition  of  q  in  (13)  by  means  of  one-sided  differences 
might  be  expected  to  limit  the  accuracy  of  this  scheme,  the  continuity  condition 
[q]  =  0  can  be  seen  to  lead  to 

0  =  [q]i=  -K"i(Ti+i -T(ei)/Ax) 

+  K+i(T(ei)-Ti_]  )/Ax  (15) 

which  indicates  that  these  flux  values  will  also  provide  second  order  accuracy  at 
interval  endpoints  which  lie  interior  of  the  solution  domain  whenever  k  is  smooth 
(in  this  case,  a  somewhat  remarkable  fact,  proved  in  [3],  is  that  the  one-sided  flux 
values  as  defined  by  (13)  are  also  2nd  order  accurate  at  the  endpoints  of  the 
underlying  domain). 

If  we  now  let 
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(a) 


(b) 


F*i=  -(qi+i-qi-i)/(2Ax)  (16) 

=  (K-i+i(Ti+i-T(ei))) 

-K+i_l(T(ei)-Ti_i  ))/(2  Ax  2) 

i  =  1,  3,...,  2M-1, 

F*i  =  -  [qli 

=  (  K+i  (T(ei+i)  -Ti)  -  K-i(Ti  -  T(ei_i))  /  Ax 

i  =  2.  4,...,  2M-2, 


then  the  use  of  the  discretization  equations  (7b),  (12),  (13),  and  (14)  leads  to  the 
following  scheme: 

(a)  ^  (g.n+l/2_  g.n)/^  ^ 

=  0. 

(b)  p*.^ri+ll2  ^  (g.n+1  _  g.n+1/2)/^^ 

=  0. 

in  which  the  variable  q  has  been  eliminated.  The  solution  of  the  tri-diagonal 
system  (17a)  yields  the  values  (Tq,  ej,  T2, 63,...,  e2j4.|, 

equations  F*i  occuring  in  (b)  also  occur  in  (a)  and  are  thus  redundant  in 

the  calculation  of  the  enthalpy  values  in  (17b).  The  nonlinearity  of  this 
system  arises  directly  through  the  terms  T(ej)  in  Fj. 

Following  [3],  we  call  (17)  the  potential  form  of  the  compact  scheme 
and  it  is  evident  from  (16)  that  F*  gives  rise  to  a  symmetric  operator.  A  proof 
(omitted  here)  that  F*  is  also  positive  can  be  given  on  the  basis  of  energy 
arguments  related  to  (4),  following  methods  described  in  [3].  The  potential  form 
thus  allows  us  to  treat  the  problem  solely  in  terms  of  the  variables  e  and  T. 


i  =  1,  3,...,  2M  -  1 
i  =  2,  4, . 2M-2 

i  =  1,  3,...,  2M  -  1 
i  =  2,  4 . .  2M-2 


(17) 
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4  A  computational  method  for  the  one-phase  Stefan  problem 

In  the  following  we  shall  assume  units  such  that  k  =  Y  =  L  =1.  Then 
7(e)  in  (7b)  has  the  specific  form 

7(e)  =  1,  0<e.  (18) 

=  0,  -1  <  e  <0 

while  (16)  reduces  to 

F*i  =  -  (Ti+i  -  2T(ei)  -  Ti_i)  /  (2  Ax  2)  (19a) 

or  else 

F=^i  =  -  (T(ei+i)  -  2Ti  -  T(ei_i))  /  Ax  (19b) 

depending  upon  the  index  i,  and  in  which 
T(e)  =  7(e)  e. 

We  see  that  (17a)  leads  to  a  nonlinear  system  of  tridiagonal  difference  equations 
involving  (T,  e)  and  to  which  are  to  be  added  boundary  conditions  involving 
T.n+1/2  and  qjn+1/2  j  ^  q,  2M. 

For  fixed  n,  if  U  indicates  the  unknowns  (T,  e),  the  symbolic  form  of 
this  algebraic  system  is 

A(U)  U  =  f  "  (20) 

which  we  propose  solving  by  the  iteration  scheme:  with  u(ll^  =  let 
solve  A(U(1^))  =  f  k=0,l,2,....  If  convergence  occurs,  then 

U  =  (T*^  +1/2  ^  ^n  +1/2^  Qyj.  previous  discussion  indicates  tliat  A(U)  is  both 
symmetric  and  positive-definite  and  thus  satisfies  conditions  which  insure  the 
solvability  of  (20)  by  an  iteration  scheme  of  this  type. 
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5  Numerical  experiments. 

Figures  1-4  illustrate  numerical  features  of  the  solution  method  just 
described.  In  the  general  study  a  parameter  o  was  used  to  relate  the  time  step  to 
the  mesh  by  means  of  At  =  o  Ax.  The  use  of  a  least  squares  fit  to  the  numerical 
position  of  the  front  verified  that  the  front  position  was  determined  with  0(At) 
accuracy.  The  average  number  of  iterations  required  to  solve  (22)  at  each  time 
step  was  in  the  range  2-4. 

The  front  position  was  determined  as  follows:  Call  a  cell  (a)  partially 
melted  if  either,  but  not  both,  the  left  or  right  endpoint  values  of  T  are  positive  , 
or  (b)  completely  melted  if  the  centerpoint  value  T(e)  >  0.  Otherwise,  call  the 
cell  unmelted.  Let 

5  =  Ax  (21) 

-  Ax/2 
=  0 

according  as  the  cell  is  completely  melted,  partially  melted  or  unmelted.  Then 
Fr(8)  =  X  5j  provides  a  distance  function  which  can  be  used  to  determine  the 

position  of  the  front.  (From  (19b)  we  observe  that  an  endpoint  value  T  >  0  iff  at 
least  one  of  the  contiguous  cells  having  this  endpoint  in  common  is  a  completely 
melted  cell). 

1.  A  front  moving  with  constant  speed: 

If,  initially,  e(x,0)  =  -I,  then  e(x,t)  =  exp(t -  x)  -  1  is  an  analytic 
solution  of  a  one-phase  Stefan  problem  for  a  slab  whose  front  moves  with  the 
constant  speed  s'(t)  =  1  when  heated  from  the  left  boundary.  This  example  thus 
allows  us  to  compare  the  numerical  results  with  known  analytic  results. 


Figures  la-lc  illustrate  results  for  the  values:  M  =  10,  a  =  1  so  that  Ax 
=  At  =  0.1.  The  choice  of  the  value  a  =  1  was  made  in  order  to  examine  a 
situation  in  which  the  effects  of  a  dispersion  error  arising  when  the  front  position 
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does  not  lie  on  grid  points  is  absent.  The  relative  errors  in  the  temperature  T  and 
the  flux  q,  using  L2  error  norms,  was  found  to  be  6.75  XlO“^  and  8.07  xl0~^ 

respectively  while  the  relative  error  of  the  front  position  was  5.09  Xl0“^^.  The 
average  number  of  iterations  per  time  step  was  2.75.  The  cause  of  the  oscillation 
in  the  q-error  in  Fig.lc  will  be  examined  shortly. 

For  M  =  20,  <j  =  1,  so  that  Ax  =  At  =  0.05.,  the  relative  errors  in  the 

temperature  T  and  the  flux  q  were  1.79  Xl0~^  and  3.92  Xl0“^  ,  respectively, 
while  the  relative  error  of  the  front  position  was  negligible.  The  average  number 
of  iterations  per  time  step  in  this  case  was  2.875. 

t  Numerical  Front 

It 


0.  8 


0.  6 


0.  4 
0.2 


— ^ - 1 - 1 - 1 - 1  X 

0.2  0.4  0.6  0.8  1 


Figure  la.  Plot  of  the  numerically  calculated  front  position  (a  =  1), 
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Figure  lb.  Plot  of  the  corresponding  numerical  temperature  T 
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E  Numerical  Solution  Error  E 
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E  Flux  Error  E 


0.  015- 
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• 
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0.  005- 
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i 

-0.01- 

• 

•  • 

Figure  I  c.  Plots  of  the  numerical  errors  of  T  and  q  (observe  the 
oscillation  in  the  q  error). 


Finally,  Figure  Id  illustrates  the  front  positions  when  Ax  =  0.1,  At  = 
0.05,. which  corresponds  to  M  =  10,  a  =  0.5.  The  relative  error  norm  of  the 

front  was  0.0542,  the  corresponding  error  norms  for  T  and  q  were  8.94  xl0“^ 

and  4.16  xl0“^,  respectively,  and  the  average  numer  of  iterations  per  time  step 
was  2,31.  As  a  result,  we  can  expect  that  the  dispersion  error  arising  from  the 
use  of  a  fixed  grid  when  a  front  with  variable  speed  occurs  will  reduce  the 
overall  accuracy.  However,  for  this  and  the  problems  considered  below  we 
found  the  accuracy  of  the  front  position  to  be  at  least  0(At). 
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2.  The  slab,  cylinder  and  sphere  heated  with  a  constant  flux  q(1.0  =  -1. 

The  preceding  results  suggest  that  the  parameter  o  may  influence  the 
accuracy  when  the  front  speed  is  variable.  This  situation  arises  in  the  following 
examples.  \ 

Figures  2-4  illustrate  numerical  results  obtained  for  various  geometries 
when  a  constant  flux  of  heat  is  added  through  at  boundary  x  =  1,  q(l,t)  =  -1, 
while  q(0,t)  =  0.  The  initial  position  of  the  front  is  s(0)  =  1.  In  all  cases,  a  =  1. 
The  average  number  of  iterations  per  time  step  was:  slab  =  2.68,  cylinder  = 
3.46,  sphere  =  3.33. 

Slab 


Figure  2a  shows  a  least  squares  fit  of  the  form 
s(t)  =  0.988  -0.883  t  +  0.118t2 

to  the  numerically  determined  front.  The  solution  plots  in  Figures  2b-2c  are  at  a 
time  slightly  after  the  time  of  the  final  front  position  indicated  in  Figure  2a.. 
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X  Numerical  Front 


Figure  2a.  A  least  squares  fit  to  the  numerical  front  for  a  slab 
geometry. 
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Figure  2b.  The  numerical  temperature  T  for  a  slab  geometry. 
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Cylinder 

The  front  positions  resulting  in  the  complete  melting  of  a  cylinder 
indicated  in  Figure  3a  has  been  fitted  with  the  cubic  polynomial 
s(t)  =  1.02  -  1.46  t  +  1.78  t  2-3.38  t^. 

The  solution  plots  in  Figures  3b  shows  the  numerical  temperature  T  prior  to  the 
time  complete  melting.has  occurred. 


X  Numerical  Front 


Figure  3a.  A  least  squares  fit  to  the  numerical  front  for  a  cylinder. 
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Figure  3b.  The  numerical  temperature  for  a  cylinder  at  a  time  shortly 
before  complete  melting. 


Sphere 


The  cubic 

s(t)  =  1.02  -  1.36  t  +  1.31  t  2-2.80  t^. 

fits  the  data  indicated  in  Figure  4a.  The  slight  irregularity  of  the  solution  which 
is  illustrated  in  Figure  4b  appears  to  be  due  to  a  reflection  arising  from  the  free 
boundary  at  x  =  s(t).  The  numerical  flux  (Fig.  4c)  appears  to  be  a  sensitive 
indicator  of  this  irregularity. 
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0^1  0^2  0^3  0^4  oTb  0l6  ^ 


Figure  4a.  A  least  squares  fit  to  the  numerical  front  for  a  sphere. 
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Figure  4b.  The  numerical  temperature  for  a  sphere  at  a  time  shortly 
before  complete  melting.  Note  the  evidence  of  an  inaccuracy  behind  the 
front. 
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Figure  4c.  The  numerical  flux  corresponding  to  Fig.  4b. 


The  parameter  o 

As  suggested  above,  the  choice  of  the  parameter  a  =  At/Ax  can  be 
expected  to  affect  the  dispersive  error  arising  whenever  the  speed  of  the  front 
differs  from  1/a.  The  following  figures  help  illustrate  this  fact  by  comparing  the 
results  for  the  sphere  shown  in  Fig.  4  (in  which  a  =  1)  with  similar  results  when 
a  =  0.2.  Figs.  5b-5c  indicate  that  this  choice  of  a  results  in  a  substantially 
smoother  approximation  and  leads  us  to  conjecture  that  a  CFL-type  condition  of 
the  form 

a  II  s'  II  <  1  (22) 

is  required  to  insure  smoothness,  rather  than  stability,  for  this  method. 
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X  Nvtmerical  Front 


Figure  5  a.  Numerical  front  for  the  sphere  when  o  =  0.2 
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Figure  5b.  Numerical  temperature  for  the  sphere  when  o  =  0.2. 
showing  a  partial  elimination  of  the  irregularity  behind  the  front 
(compare  Figure  4b). 
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Figure  5c.  Partial  smoothing  of  the  numerical  flux  when  o  =  0.2 
(compare  to  Fig.  4c). 
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Concluding  Remarks: 

These  results  indicate  that  the  ID,  one-phase  Stefan  problem  can  be 
treated  by  the  enthalpy  method  using  an  unconditionally  stable  numerical  method 
which  requires  the  iterative  solution  of  a  symmetric  (positive-definite)  tridiagonal 
system  of  equations  (21)  at  each  time  step.  Because  the  number  of  iterations  per 
time  step  is  relatively  small  (<  4  for  each  of  the  cases  treated),  this  provides  a 
computationally  efficient  method.  Also,  because  this  same  basic  structure 
underlies  the  use  of  ADI  solution  techniques  it  suggests  that  multidimensional 
problems  can  be  treated  in  a  similar  manner. 

Further  studies  are  required  to  understand  the  mechanisms  which  give 
rise  to  numerical  irregularities  when  the  front  moves  relative  to  a  fixed  mesh  and 
which  appear  to  be  controlled  by  the  parameter  a  =  At/Ax.  Other  aspects  of  this 
problem  can  be  expected  to  appear  in  higher  dimensions. 

Finally,  in  view  of  the  fact  that  discrete  "energy"  estimates  similar  to 
those  available  for  smooth  solutions  as  described  by  eq.(5)  can  be  given  for  tlie 
scheme  described  here  (following  developments  given  in  Rose  [4]),  the  possibility 
exists  of  developing  a  more  complete  theoretical  analysis  of  the  enthalpy  method 
based  on  this  observation.  In  particular,  this  would  prove  that  the  basic  operator 
arising  from  the  scheme  is,  indeed,  positive  definite,  a  fact  that  was  not  shown  in 
this  paper. 
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Abstract.  We  consider  both  the  compact  finite  volume  and  finite  difference  space 
discretizations  of  the  Stefan  problem.  The  resulting  algebraic  systems  are  solved  by  nonlinear 
versions  of  ADI  and  SOR.  Both  algorithms  contain  significant  parallelism  which  is  demonstrated 
on  two  vector/multiprocessing  computers,  the  Alliant  FX/40  and  the  Cray  Y-MP.  Numerical 
experiments  indicate  that  the  compact  discretization  and  ADI  give  the  best  accuracy  with  the 
minimum  computational  cost. 
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A  COMPACT  FINITE  VOLUME  SCHEME  FOR  2-D  STEFAN 
PROBLEMS  AND  VECTOR/MULTIPROCESSOR  COMPUTERS. 


1.  Introduction.  This  report  outlines  some  new  numerical  methods  for  the  solution 
of  multiple  space  dimension  Stefan  problem.  This  will  include  both  SOR  and  ADI  algorithms  for 
the  nonlinear  algebraic  systems  which  result  from  both  the  FDM  and  compact  space  discretizations 
methods.  Analysis  of  convergence  will  be  discused. 

The  ADI  algorithms  have  been  introduced  because  they  appear,  in  many  examples,  to  be 
more  effective  than  some  of  the  traditional  algorithms  related  to  the  Stefan  problem.  Moreover, 
these  ADI  algorithms,  as  in  the  linear  parabolic  problems,  have  large  number  of  independent 
tridiagonal  solvers.  Therefore  they  can  be  effectively  implemented  on  vector/multiprocessing 
computers.  This  is  illustrated  for  Alliant  FX/40  and  the  Cray  Y-MP. 

In  this  paper,  we  will  use  the  enthalpy  formulation  of  the  Stefan  problem.  (See  Rose  [3], 
Elliot  [2]  and  White  [5]). 

Et  -  Ap(E)  =  f  on  D  X  (0,T) 

E(x,0)  =  given  on  D 

(1)  P(E(x,t))  =  given  on  dD  X  (0,T), 

where  P(E)  =  Kirchoff  transformed  temperature 

P(E)  =  u  is  the  “inverse”  of  the  enthalpy  function  E  =  H(u). 
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Both  (8.1)  and  (8.2)  can  be  solved  as  in  (5).  For  example,  consider  (8.1)  with  k  +  1/2 
suppressed  and  m  equals  the  iteration  index. 

Consider 


(9)  (aI  +  H(E))F.  =  d 

E*"  + 1  =  (a  I  +  H(E"'))  ‘  d 

(10)  =(I+-l-H(E"'))*a-*d 

a 


Proposition  1: 

Let  A,  A(E),  H(E),  V(E)  be  as  above.  If  A  is  an  M-matrix  and  P(E)  is  as  defined  above,  then 

A(E),  H(E),  V(E)  are  non-singular  for  small  At.  Moreover,  there  exists  Oo  >  0  such  that  for  a  >  Oq 

(10)  converges  to  the  unique  solution  of  (9). 

The  proof  follows  that  given  in  White  [6]  for  (4)  and  (5). 

In  practice  a  is  used  as  an  acceleration  parameter  and  usually  only  3-5  iterations  are 
required  to  solve  (8.1)  or  (8.2).  The  schemes  in  (8.1)  and  (8.2)  do  not  solve  the  implicit  time  step: 

(11)  +  +  =  +  ‘ 

At 

For  each  fixed  problem  in  (1 1)  there  are  a  number  of  ADI  schemes  with  multiple  H  and  V 
sweeps.  We  focus  on  only  one  in  the  next  section. 
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3.  ADI:  Multiple  H  and  V  sweeps  per  time  step.  Consider  the  simple 
nonlinear  system  (4)  with 


A(E)  =  H(E)  +  V(E),  or 
(12)  (H(E)  +  V(E))E=d 

After  linearization  the  multiple  H  and  V  sweep  ADI  algorithm  is 

(13.1)  (a  I  +  H(E'^))  E'"  +  1/2  =  +  (ct  I .  v(E'^))  E*" 

(13.2)  (a  I  +  V(E‘^))  E"*  ^  1  =  d  +  (a  I  -  H(E"’'>))  E*"  ^ 

One  can  solve  (13.1)  and  insert  it  into  (13.2)  to  obtain 

(14)  E'"  + 1  =  5(E'^)  +  JKE"^)  E"™ 

If  E"*  1  converges,  defuie  E"'®  ^  i  to  be  the  limit 

Etno  +  1  =  5(E‘"®)  +  HfE""®)  E*^  +  1 
provided  the  spectral  radius  of  HfE"*®)  is  less  than  1.  This  is 

(15)  E"’°  ^  1  =  a  *  H(E">®))‘‘  5(E"'®) 

Line  (15)  forms  the  outer  iteration  of  our  algorithm  to  solve  (2). 

Proposition  2: 

Consider  (12)  and  the  algorithms  (13)  -  (15).  Let  the  same  assumptions  hold  as  in  proposition 
1.  Then  there  exists  a  Oq  >  0  such  that  if  a  >  Oq,  then  the  inner  iterations  (14)  converges  to  E*"®  ^  i 

in  (15),  and  E*”®  ^  i  in  (15)  converges  to  thr  solution  in  (12). 

The  proof  requires  the  spectral  radius  of  H(E)  to  be  uniformly  bounded  below  1.0.  This 
step  makes  use  of  monotonicity  properties  (nonnegative  matrices),  and  it  contrasts  with  the  linear 
case  where  symmetric  positive  definite  matrices  are  used. 

The  above  schemes  is  not  the  only  way  to  do  multiple  H  and  V  sweeps,  but  it  does  allow 
one  to  give  a  convergence  analysis.  The  best  version  of  multiple  H  and  V  sweeps  is  not  yet  clear. 
The  numerical  experiments  that  we  have  done  indicate  that  the  Peaceman-Ranchford  method  (one 
H  and  V  sweep  per  time  step)  is  adequate. 


8 


4.  Numerical  Experiments.  In  this  section  we  consider  the  numerical  solution  of 


Et  -  AP(E)  =  0,  where  the  domain  D  =  (0,1)  x  (0,T),  T  = 


The  space  variables  are  discretized  in  two  ways,  the  compact  finite  volume  (CRT)  (See 
Rose  [4]),  and  the  traditional  five-point  finite  difference  method  (FDM).  The  resulting  algebraic 
problems  are  solved  by  ADI  (one  H  and  V  sweep  per  time  step)  and  by  SOR.  The  computations 
were  done  on  two  vector/multiprocessing  computers:  the  AUiant  FX/40  and  the  Cray  Y-MP.  The 
Alliant  FX/40  has  two  processors  each  with  a  vector  pipeline.  The  Cray  Y-MP  was  used  as  a 
single  processor  with  a  more  sofisticated  vector  pipeline  and  much  shorter  clock  cycle  time  (4  nsec 
as  compared  to  170  nsec). 

The  CPT  space  discretization  has  a  great  deal  of  parallelism  wich  can  be  executed  on  either 
vector  pipelines  or  multiprocessing  architectures.  The  unknowns  at  the  center  of  the  cells  maybe 
grouped  together,  and  the  unknowns  at  the  boundary  of  the  cells  grouped  together.  The  resulting 
coefficient  matrix  has  the  form  of  a  two  coloring  scheme.  Hence  one  can  execute  the  SOR  scheme 
for  the  CPT  discretization  on  vector  pipelines.  Also  the  CPT  space  discretization  when  solved  by 
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ADI  methods  will  have  a  large  number  of  independent  tridiagonal  solvers.  Thus,  the  vector 
version  of  the  tridiagonal  algorithm  may  be  used. 

In  the  case  of  FDM  space  discretization,  the  traditional  red-black  ordering  of  nodes  allows 
vectorization  of  the  SOR  algorithm.  Also,  domain-decomposition  methods  have  been  used  on 
SOR  so  that  the  two  processors  of  the  Alliant  FX/40  can  be  effectively  used.  Of  course,  we  can 
use  the  vector  version  of  the  tridiagonal  algorithm  to  solve  the  FDM  when  ADI  is  used. 


V.  CeUs 

Computers 

2 

16 

32^ 

2 

64 

Alliant,  serial 

2.01 

15.64 

122.28 

Alliant,  vector 

1.02 

5.98 

38.63 

Alliant,  vector,  2  processors 

0.57 

3.32 

21.07 

Cray,  serial 

0.129 

0.89 

6.52 

Cray,  vector 

0.062 

0.32 

1.76 

Table  1: 

CPT-ADI 

Cells 

Computers 

2 

16 

2 

32 

2 

64 

Alliant,  serial 

4.68 

74.44 

1155.00 

Alliant,  vector 

2.67 

30.72 

508.81 

Alliant,  vector,  2  processors 

1.41 

15.99 

270.91 

Cray,  serial 

0.41 

6.23 

96.24 

Cray,  vector 

0.11 

1.06 

11.34 

Table  2: 

CPT-SOR 
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Table  1  contains  computing  time  done  for  the  compact  discretization  and  using  one  sweep 
ADI  algorithm.  Both  computers  were  effective  vectorizers.  Table  2  contains  computing  times 

done  for  the  compact  discretizatin  with  red-black  ordering  SOR  algorithm  where  O)  =  1.2. 


Cells 

Computers 

2 

16 

2 

32 

2 

64 

Alliant,  serial 

0.82 

12.93 

186.72 

Alliant,  vector 

0.89 

10.29 

124.97 

Alliant,  vector,  2  processors 

0.51 

5.58 

67.35 

Cray,  serial 

0.069 

1.05 

15.48 

Cray,  vector 

0.043 

0.42 

4.19 

Table  3:  FDM-SOR 


The  third  table  contains  cc  mputing  times  for  the  FDM  with  red-black  ordering  for  SOR 
where  co  =  1.2.  This  method  requires  more  iterations  to  reach  convergence  than  does  the  CPT- 
SOR  method. 

All  three  methods  give  good  accuracy.  The  CPT-ADI  seems  to  be  more  accurate  and  less 
computing  time  required  for  this  and  related  examples. 
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Abstract.  Recently  developed  Compact  Finite  Difference  scheme  (CPT)  is  applied  to 
two  dimensional  diffusion  equations.  The  relative  merits  of  CPT-ADI  are  investigated  with  other 
computational  schemes  such  as  finite  difference  method  -  ADI  (FDM-ADI)  and  FDM-SOR.  The 
numerical  results  obtained  from  these  three  approaches  are  compared  to  known  analytical 
solutions.  The  primary  intere.st  of  this  study  lies  on  vectorization  and  parallel  processing. 
According  to  our  results  shown  in  tables  !,  2  and  3,  CPT-ADI  is  found  to  be  superior  scheme  with 
regards  to  accuracy,  than  both  FDM-ADI  and  FDM-SOR.  It  is  al.so  fastest  algorithm  than  both 
FDM-ADI  and  FDM-SOR  as  it  is  evident  from  CPU  times. 


*  This  research  was  sponsored  by  the  Air  Force  Office  of  Scientific  Reseiu'ch,  Bolling  Air 
Force  Base,  D.  C.  Contract  No.  F49620-89-C-001()  and  Anny  Research  Office,  Research  Triangle 
Park,  N.  C.  Contract  No.  DAAL03-90-G-0126. 


A  COMPARATIVE  STUDIES  OF  COMPACT  FINITE  VOLUME 
METHODS  FOR  THE  2-D  DIFFUSION  EQUATION 
WITH  FINITE  DIFFERENCE  ADI  AND  SOR. 


1.  Introduction.  We  shall  describe  briefly  the  compact  finite  difference  method 
(CFT)  for  one  dimmensional  steady  state  problem.  The  extension  to  2-D  problem  may  be  easily 
done.  The  underlying  physics  behind  this  approach  lies  in  solving  the  differential  equation  in 
isolation  from  its  neighboring  subintervals  (i.e.  compactiy).  Then,  extend  the  solution  in  the  large 
by  means  of  continuity  conditions  for  the  flux  and  temperature  accross  the  boundaries  of  the 
contiguous  subintervals. 

We  consider  here  one  dimensional  steady  diffusion  problem: 


(1) 

D  F. 

u'  =  f 

(2) 

(j)"  =  u 

(3) 

B.  C. 

4>  =  g 

for  x  €  I 

and  I=[I.,  14 

Divide  the  interval  I  into  m  nonoverlaping  subintervals: 

Ij  =  X  I  X-  1<  X  <  Xj^l), 

I  2  21 


Axj 


(Axj  =  2  hj) 


1  ^  1 

with  center  points  xj,  J  ~  2  ’  2 ’  ’  o  interior  endpoints  xj,  j  =  1,2, ..., 


adopt  the  finite  difference  notations  Axj  =  x^  ^  -  x^  k  hj  =  ~  and  u(i)  =  Ui. 

2  2 

Ic  =  1^,  2 . ^  '  2 1  points) 

le  =  {1,  2,  ...,  M  -  1)  (interior  endpoints) 


M  -  1.  We  shall 
We  also  denote: 


K 


Ax 


Ax 


2  2  2 

Fig.  1. 

The  discrete  equations  for  (1)  and  (2)  can  be  written  as,  for  j  e  Ic, 


(4) 

2  2 

(5  a) 

2  2 

(5  b) 

and  0;  -  (j); .  1  =  h;  U-  1 

2  ■'‘2 

Next,  it  is  required  that  both  u  and  <j)  be  continuous  across  every  endpoint  common  to  two 
intervals: 

[u]i  =  [4  =  0  for  i  e  Ig. 

Using  this  continuity  condition  [u],  =  0  in  terms  of  (|>  we  have 


4>.-l 

2 


2 


_ _  2 


h,.l 


hi+i 


and  in  the  case  of  equal  intervals,  which  reduces  to 

(j)i  + 1  +  (l)i .  1 


i  €  Ic  (endpoints) 
4 


(6) 


Now  from  (4)  and  (5)  we  get. 


(t)j  +  1  -  2  ({)j  +  ([ij  .1  =  2  hj  fj,  i  G  Ic 


and  from  (6) 
(7  b) 


(l)j  +1  -  2  (l)j  +  ([);  .1  =  0,  i  €  le 
2  2 


The  equations  (7)  and  the  boundary  conditions  lead  to  determined  system  of  algebraic 
equations  for  the  values  of  (J).  These  equations  lead  to  a  tridiagonal  system  of  equations  which  is, 
therefore,  solved  by  Thomas  algorithm. 

The  extension  to  two  dimensional  scheme  may  be  easily  done.  Two  dimensional  stencil  is 
shown  in  Fig.  2. 


Ui,j+  1 


".  +  i.j+J- 

7 


Ui  +  l.j 


Consider  the  two  dimensional  diffusion  equation  Ut  =  Uxx  +  Uyy-  The  compact  ADI 
method  is  given  by  lines  (8)  -  (12): 


n  +  J- 


(8  a) 

(8  b) 


u  2-U?.j 

- =  F  2  +  G"  : 

X  i.  j  • 


n  +  - 


if?  t  ■  I*  2  1 

- LJ_=f  I  + 

t  i.  j  '•  J 

where  F;,  j  and  G;,  j  are  standard  finite  difference  expressions  for  u^x  and  Uyy  respectively  and 


X  =  Now 
2 


n  +  J- 


(9  a) 

F  2  = 
i.j 

(9  b) 

c^.j  = 

n  +  J-  n  +-1- 

.2: 


2  Ax 


2  Ay 

Notice  that  we  modify  the  indices  from  the  previous  discusion.  From  Fourier  law  we  have 
q  =  -KUx,  therefore 


(10  a) 

(10  b) 
(10  c) 

(10  d) 


q  2  = 

"  +1.J 


h-i.j 


n  +-!- 

n  +-!- 

+  i.j 

U  2 

-  U  2 

.  1+  I.j 

'.J 

Ax 

n  +1- 

n  +  J- 

-Ki-l.j 

U  2  - 

.  '.j 

U  2 

i-l.jJ 

Ax 

ik.j^ 

Ay 

l[<j- 

<j-i] 

Ay 

Assuming  K  =  1  and  Ax  =  Ay,  we  have  from  (8)  with  the  help  of  (9)  and  (10) 
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(11) 


n  +  — 


n  +—  n  +  — 


u.  2 .  +  (a  +  2)  u.  .  2  -  u.  2  =  u?  J .  1  +  (a  -  2)  j  -  u?  j  1 

I  -  i,j  1,  j  j  +  i.j  *' j  *  ’  j  ’j 


i,  j  =  1,  3, 2  M  -  1,  a  =  ^  ,  Ax  =  Ay 


At 


(12) 


1+1  „n  +  1  _ 

J  '“i.j+l 

i,  j  =  1,  3, 2  M  -  1 


n  +  1 


,  ,  .  ...  n  +  J-  n  +  J- 

-u^  j  .*1  +  (a  +  2)  j  -  u?  j  +  1  =  u.  2 .  +  (a  -  2)  u.  .  2  -  u.  ^  2  ^ 


and  the  continuity  conditions  are 


(11a) 

n+1 

n+1 

-U  2 

+  2  U  2 

-  U  2 

i  -  l.j 

‘.J 

1  +  l.j 

i,j  =  2,  4, ... 

,2M-2 

(12  a) 

+  2  <  j  J 

„n  +  1 

*  “i.j  +  1 

i,  j  =  2,  4, ... 

,  2  M-2 

The  two  algebraic  systems  (11),  (11a),  and  (12),  (12a)  can  be  solved  by  the  tridiagonal 
algorithm. 


7 


2.  Numerical  experiments.  One  of  the  primary  objectives  of  this  project  is  to 
vectorization  of  1-D  and  2-D  diffusion  problems.  We  vectorize  the  2-D  heat  equation  for  the  Cray 
Y-MP  using  red-black  SOR  algorithm. 

We  consider  three  analytical  solutions  for  the  following  diffusion  equation  and  compare 
them  with  respective  numerical  solutions.  The  results  are  recorded  in  tables  1,2  and  3. 

The  two  dimensional  diffusion  equation  is  Ut  =  Uxx  +  Uyy  +  f(x,y,t),  0<x<  l,0<y<  1, 

t  >  0.  The  following  analytical  solutions  are  considered: 

(a)  u(x,  y,  t)  =  exp{t  -  x  cos0  -  y  sin0 },  6  being  any  angle,  however,  0  =  0, 7t/2  reduce 

the  problem  into  an  one  dimensional  one.  We  use  f  =  0  and  0  =  7i/4. 

(b)  u(x,  y,  t)  =  2  t  +  x2  cos^G  +  y2  sin^G,  G  being  any  angle,  f=  0. 

(c)  u(x.  y,  t)  =  X  (1  -  x)  y  (1  -  y)e-‘,  and  f(x,  y,  t)  =  (x  (1  -  x)  (2  -  y  (1-y))  +  2  y  (l-y))e-‘. 

All  the  problems  mentioned  above  are  associated  with  Dirichlet  boundary  conditions. 

The  CPT-ADI  is  compared  with  FDM-ADI  and  FDM-SOR.  The  functional  errors  in  each 

of  the  three  methods  are  of  Ofh^)  as  expected  and  derivative  errors  are  of  0(h).  The  results  are 

shown  in  table  1,  table  2  and  table  3. 
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FDM-ADI  At  =  (Ax)^ 


No.  of 

hx  =  hy  cells 

No.  of 

iterations 

Max 

function 

error 

Max 

derivative 

error 

8 

32 

0.2834 

2.07583E-2 

1.055968E-1 

16 

128 

4.45 

6.48407E-3 

5.343229E-2 

32 

512 

69.65 

1.826602E-3 

2.889879E-2 

64 

2048 

1105.48 

8.967522E-4 

1.5545622E-2 

CPT-ADI  At  =  Ax 

No.  of 

No.  of 

CPU 

Max 

Max 

time 

function 

derivative 

hx  =  hy  cells 

iterations 

(sec) 

error 

error 

8 

8 

0.15 

9.359392E-5 

1.9731343E-3 

16 

16 

1.0 

2.366787E-5 

6.054254E-4 

32 

32 

7.567 

5.9340348E-6 

2.120109E-4 

64 

64 

58.8 

1.4842674E-6 

8.433524E-5 

FDM-SOR  At  = 

Ax 

No.  of 

No.  of 

CPU 

Max 

Max 

time 

function 

derivative 

hx  =  hy  cells 

iterations 

(sec) 

error 

error 

8 

8  X  7  =  56 

0.134 

1.20095E-2 

9.07367E-2 

0)=  1.55 

16 

3.4 

3.21928E-3 

5.0422E-2 

(0=  1.7 

32 

49.85 

8.3834E-4 

2.16456E-2 

(0=  1.7 

64 

S4  X  92  =  588^ 

1122.8 

2.3626E-4 

1.1279E-2 

(0=  1.7 

Table  1.  Problem  (a)  u(x,  y,  t)  =  exp(t  -  x  cos0  -  y  sin0),  6  =  7t/4 
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FDM-ADI  At  =  (Ax)^ 


No.  of 

hx  =  hy  cells 

No.  of 

iterations 

CPU 

time 

(sec) 

Max 

function 

error 

Max 

derivative 

error 

8 

32 

0.2834 

2.07583E-2 

1.055968E-1 

16 

128 

4.45 

6.48407E-3 

5.343229E-2 

32 

512 

69.65 

1.826602E-3 

2.889879E-2 

64 

2048 

1105.48 

8.967522E-4 

1.5545622E-2 

CPT-ADI  At  =  Ax 


No.  of 

hx  =  hy  cells 

No.  of 

iterations 

CPU 

time 

(sec) 

Max 

function 

error 

Max 

derivative 

error 

8 

8 

0.15 

9.359392E-5 

1.973 1343E-3 

16 

16 

1.0 

2.366787E-5 

6.054254E-4 

32 

32 

7.567 

5.9340348E-6 

2.120109E-4 

64 

64 

58.8 

1.4842674E-6 

8.433524E-5 

FDM-SOR  At  = 

Ax 

No.  of 

hx  =  hy  cells 

No.  of 

iterations 

CPU 

time 

(sec) 

Max 

function 

error 

Max 

derivative 

error 

8 

8  X  7  =  56 

(0=  1.55 

0.134 

1.20095E-2 

9.07367E-2 

16 

3.4 

3.21928E-3 

5.0422E-2 

(0=  1.7 

32 

49.85 

8.3834E-4 

2.16456E-2 

co=  1.7 

64 

S4  X  92  =  588^ 
(0=  1.7 

1122.8 

2.3626E-4 

i.l279E-2 

Table  2.  Pi  oblem  (b)  u(x,  y,  t)  =  2  t  +  x^  cos^G  +  y2  sin^9,  0  =  ji/4 
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FDM-ADI  At  =  (Ax)^ 


No.  of 

hx  =  hy  cells 

No.  of 

iterations 

Max 

function 

error 

Max 

derivative 

error 

8 

32 

0.4167 

1.3274E-4 

1.202978E-2 

16 

128 

6.7667 

5.55037E-5 

6.083 185E-3 

32 

512 

110.28 

2.80686E-5 

3.063236E-3 

64 

2048 

1786.05 

1.39326E-5 

1.538248E-3 

CPT-ADI  At  =  Ax 

No.  of 

No.  of 

CPU 

Max 

Max 

time 

function 

derivative 

hx  =  hy  cells 

iterations 

(sec) 

error 

error 

8 

8 

0.1834 

6.8486  lE-4 

2.73903E-3 

16 

16 

1.3167 

3.521 155E-4 

1.408662E-3 

32 

32 

10.267 

1.780065E-4 

7.121056E-4 

64 

64 

80.834 

8.942948E-5 

3.57749E-4 

FDM-SOR  At  = 

Ax 

No.  of 

No.  of 

CPU 

Max 

Max 

j 

time 

function 

derivative 

hx  =  hy  cells 

iterations 

(sec) 

error 

error 

8 

8x1=8 

5.0E-2 

2.802506E-3 

1.842519E-2 

(0=1.7 

16 

16x2  =  32 

1.1334 

1.089026E-3 

1.286985E-2 

(0=1.9 

32 

32x5  =  160 

12.784 

6.681835E-4 

4.63165E-3 

(0=  1.85 

64 

54  X  16  =  1024 

319.58 

8.712826E-5 

1.700953E-3 

(0=  1.9 

Table  3.  Problem  (c)  u(x,  y,  t)  =  x  (1  -  x)  y  (1  -  y)  e  ‘ 
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3.  Concluding  remarks.  The  x,  y  domains  are  devided  into  8,  16,  32  and  64  cells 
and  in  each  case,  the  starting  time  is  taken  to  be  0  and  the  final  time  is  1 .  The  relation  between  time 
step  and  space  step  varies  among  FDM-ADI,  FDM-SOR  and  CPT-ADI  schemes.  We  assume  all 

the  material  constants  are  equal  to  unity  and  Ax  =  Ay.  In  case  of  FDM-SOR  and  CPT-ADI,  Ax  = 

Ay  =  At  and  however,  in  case  of  FDM-ADI,  At  =  (Ax)  .  This  is  a  disadvantage  for  FDM-ADI.  For 

example,  when  Dx  =  1.5625E-2,  Dt  =  4.883E-4  and  it  would  take  2048  iterations  to  reach  time  1. 
However,  CPT-ADI  and  FDM-SOR  are  free  from  this  difficulty.  But  SOR  also  has  a  different 
disadvantage.  Each  time  step,  FDM-SOR  takes  large  number  of  iterations  to  converge  to  a 
preassigncd  level.  With  64  cells  CPT-ADI  requires  only  64  iterations  to  reduce  the  error  level  to 
4.84303E-7,  where  as  FDM-ADI  requires  2048  iterations  to  educe  the  error  level  to  4.081723E-4 
and  FDM-SOR  requires  5,952  iterations  to  reduce  the  error  level  to  7.77 14E-6. 

Perhaps  the  most  convincing  argument  may  be  made  in  favor  of  CPT-ADI  with  respect  to 
CPU  times.  In  this  respect,  CPT-ADI  is  the  best  With  64  cells,  CPT-ADI  requires  only  58.8  secs 
where  as  FDM-ADI  and  FDM-SOR  need  1105.48  secs  and  1122.8  secs  respectively  to  compute 
the  same  problem. 

Acknowledgments:  Thanks  to  Milton  E.  Rose  for  his  many  valuable  suggestions  for  the 
completion  of  this  paper.  Also  many  thanks  go  to  Arje  Nachman  of  AFOSR. 
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ATTACHMENT  #5 


-  BStefan  2D.7Report  -  2/14, ^91  - 


Stefan  2D 

This  parer  reports  on  the  solution  of  the  Stefan  problem 
in  2  dimensions.  The  codes  are  written  in  True  Basic,  and 
the  Enthalpy  formulation  of  the  2D  Stefan  problem  is 
approximated  by  Compact  Schemes.  The  niomerical  results 
are  compared  to  an  exponential  solution,  and  the  solution 
and  errors  are  plotted  using  Mathematica. 

Dirichlet  boundary  conditions 

The  following  results  pertain  to  an  exponential  solution 
of  the  Stefan  problem.  The  melting  front  moves  on  lines 

t  -  XCOS0  -  ysin6  =  K  with  speed  1. 

0  =  0 

The  parameters  set  in  the  code  specified  32  cells  in  all 
directions,  and  solution  in  [0,  .7]  time  interval.  The 

problem  is  reduced  to  ID. 

Import  data 
Plot  Data 

A  contour  plot  of  the  analytic  solution  would  reveal 
horizontal  lines  behind  the  front  at  position  22.4. 

Below  is  the  contour  plot  of  the  numerical  solution. 


-  BStefan  2D.7Repoil  -  2/14/91  - 


ListContourPlot [Cl,  PlotLabel -> "Temperature " , 
ContourLevel a - >2  0 ] 

30- 

25 

2  0- 

15- 

10- 

5 

5  10  15  20  25  30 

-ContourGraphics- 

A  3D  plot  of  the  numerical  solution  shows  the  behaviour 
in  the  neighborhood  of  the  front. 


-2- 


BStefan  2D.7Report  -  2/14/91  - 


Li8tPlot3D[Sl,  PlotLabel-> "Temperature",  Shading- >False 
ViewPoint->{4.000,  0.910,  1.290}] 

Temperature 


-SurfaceGraphics- 

Below  is  a  plot  of  the  positions  of  the  numerical  front 
for  each  corresponding  time  frame. 


-  BStefan  2D.7Report  -  2/14/91  - 


ListPlot3D [Fl,  PlotLabel-> "Front  tracking".  Shading- >False, 
ViewPoint-> { 0 .790,  -2.570,  1.290>] 

Front  tracking 


-SurfaceGraphics- 

A  3D  plot  of  the  numerical  error  shows  that  it  is  mostly 
confined  to  the  neighborhood  of  the  front. 


-4- 


-  BStefan  2D.7Report  -  2/14/91  - 


LlstPlot3D [El,  PlotLabel-> "Error",  Shading- >Fal8e, 
PlotRange->All,  ViewPoint-> {4 . 000,  0.910,  1.290}] 

Error 


-SurfaceGraphics- 


0  =  7i:/4 


The  parameters  set  in  the  code  specified  32  cells  in  all 
directions,  and  solution  in  [0,  .7]  time  interval. 

S  Import  data 


BPIot  Data 

A  contour  plot  of  the  analytic  solution  would  reveal 
diagonal  lines  behind  the  front.  Below  is  the  contour 
plot  of  the  numerical  solution. 
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-  BStefan  2D.7Report  -  2/14/91  - 


ListContourPlot [C2,  PlotLabel- > "Temperature " , 
Cont ou r Le ve 1 s - > 2 0 ] 


Temperature 


A  3D  plot  of  the  numerical  solution  shows  the  behaviour 
in  the  neighborhood  of  the  front. 
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-  BStefan  2D.7Repon  -  2/14/91  - 


Li8tPlot3D[F2,  PlotLabel->"Front  tracking",  Shading->False, 
ViewPoint->{0.790,  -2.570,  1.290}] 

Front  tracking 


-SurfaceGraphics- 

A  3D  plot  of  the  numerical  error  shows  that  it  is  mostly- 
confined  to  the  neighborhood  of  the  front. 
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LiStPlot3DlE2 
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■  import  doto 
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-  BStefan  2D.7Report  -  2/14/91  - 


■  Plot  Data 

Below  is  the  contour  plot  of  the  ntimerical  solution. 

ListContourPlot [C4,  PlotLabel -> "Temperature " , 
ContourLeve 1 s - >2  0 ] 


Temperature 


A  3D  plot  of  the  numerical  solution  shows  the  behaviour 
in  the  neighborhood  of  the  front. 
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-  BStefan  2D.7Report  -  2/14/91  - 


ListPlot3D[S4,  PlotLabel->"Teinperature" ,  Shading->False, 
ViewPoint->{4.000,  0.910,  1.290)] 

Temperature 


-Surf aceGraphics- 

Below  is  a  plot  of  the  positions  of  the  numerical  front 
for  each  corresponding  time  frame. 
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-  BStefan  2D.7Report  -  2/14/91  - 


ListPlot3D[F4,  PlotLabel-> "Front  tracking",  Shading->False, 
ViewPoint-> {0.790,  -2.570,  1.290}] 

Front  tracking 


-SurfaceGraphics- 

A  3D  plot  of  the  nvimerical  error  shows  that  it  is  mostly 
confined  to  the  neighborhood  of  the  front. 
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-  BStefan  2D.7Report  -  2/14/91  - 


ListPlot3D[E4,  PlotLabel->"Error",  Shading- >False, 
PlotRange->All,  ViewPoint->{4 . 000,  0.910,  1.290}] 

Error 


-SurfaceGraphics- 

Rll  FIuh  boundary  conditions 

0  =  0 


0  =  7c/4 


The  parameters  set  in  the  code  specified  32  cells  in  all 
directions,  and  solution  in  [0,  .7]  time  interval. 

■  Import  data 

■  Plot  Data 

Below  is  the  contour  plot  of  the  numerical  solution. 
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-  BStefan  2D.7Report  -  2/14/91  - 


ListContourPlot [C6,  PlotLabel -> "Temperature " , 
ContourLevel a - >20 ] 


Temperature 


-ContourGraphics- 

A  3D  plot  of  the  numerical  solution  shows  the  behaviour 
in  the  neighborhood  of  the  front. 
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BStefan  2D.7Report  -  2/14/91  - 


ListPlotBD [S6,  PlotLabel-> "Temperature",  Shading- >False 
ViewPoint->{4.000,  0.910,  1.290}] 

Temperature 


-SurfaceGraphics- 

Below  is  a  plot  of  the  positions  of  the  numerical  front 
for  each  corresponding  time  frame. 
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BStefan  2D.7Report  -  2/14/91  - 


ListPlot3D [F6,  PlotLabel-> "Front  tracking".  Shading- >False 
ViewPoint->{0.790,  -2.570,  1.290}] 

Front  tracking 


-Surf aceGraphics- 

A  3D  plot  of  the  numerical  error  shows  that  it  is  mostly- 
confined  to  the  neighborhood  of  the  front. 
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ListPlotSD [E6,  PlotLabel->''Error",  Shading- >False 
PlotRange->All,  ViewPoint->{4 . 000,  0.910,  1.290}] 

Error 


-SurfaceGraphics- 


